Season of prescribed burns and management of an early successional species affect flower density and pollinator activity in a pine savanna ecosystem

Study design and site description

This study occurred at Austin Cary Forest (ACF) (29°43′N, 83°13′W), a 2,080-acre forest in Alachua County, Florida, USA. ACF was historically used as a timber and turpentine operation before the University of Florida purchased the area in 1930; now, most of the land is used for research and teaching. ACF is managed with prescribed burns that vary in frequency and season.

The study plots were embedded in a 10-hectare early successional longleaf pine savanna. The broader study area is classified as a xeric upland pine savanna with predominantly Entisol soils. Upland pine savannas such as our study site are typically ignited every 2–5 years (Huffman, 2006). The entire area had previously been logged in 2013, burned in 2016, and planted in 2017 with longleaf pines, Pinus palustris. These restoration actions caused soil disturbance, which resulted in the dominance of yankeeweed, a perennial, ruderal species common in early successional fields. In late 2017, the area was divided into 12 plots and assigned to seasonal burn treatments, varying in the timing of prescribed burns, as described below. At the start of this study, the time since fire was identical across all plots, which had been burned in 2016 prior to initiating this experiment.

Timing of prescribed burn treatments

We defined our burn treatments using precipitation patterns to easily compare our results with other study systems (Fill, Davis & Crandall, 2019; Baruzzi et al., 2022). Based on cumulative rainfall anomalies, which indicate trends of increasing and decreasing precipitation based on long-term averages, the timing of prescribed burn aligned with three distinct seasons: (1) winter, dry season burns hereafter referred to as winter-dry season burns; (2) late spring, early-wet season burns, hereafter referred to as spring burns; and (3) summer, wet season burns, hereafter referred to as summer-wet season burns. These three burn treatments and an unburned control were replicated three times across 12 plots. Treatments were randomly assigned to plots, with the exception that if random selection resulted in a treatment occurring in two adjacent plots, another treatment was randomly selected. Each burn treatment plot was approximately one hectare in size, but the exact plot size varied to allow for the placement of adequate fire breaks between plots. Burn treatment plots were directly adjacent to one another but with a fire break in between each plot. Seasonal burn treatments were first applied in 2018 prior to data collection. Winter-dry season burns were applied on March 22, 2018, spring burns on May 25, 2018, and summer-wet season burns on August 16, 2018. Data on flowers and pollinators were collected throughout 2019. Prescribed burns were conducted again in 2020, with winter-dry season burns applied on January 29, 2020, spring burns on June 18, 2020, and summer-wet season burns on August 9, 2020, and data were collected on flowers and pollinators in that same year (Fig. 1). For all seasonal burn treatments, backing fires were lit against the wind to establish a blackline before lighting successive strip fires across each plot. If subplots did not burn, the edges were relit with a head fire to ensure most of all plots were burned. Generally, summer-wet season burns consumed the most biomass with >95% of fuel burned, followed by winter-dry and spring season burns, which ranged from 75–95% fuel consumption across plots.

Yankeeweed removal treatments

Within each fire treatment plot, two 10 m × 10 m yankeeweed treatment subplots were established, including one control subplot in which yankeeweed was not manipulated, hereafter referred to as the control plot, and one removal subplot in which yankeeweed was removed (~100% coverage reduction), hereafter referred to as the removal plot. Each subplot was located at least 10 m from plot borders. Within yankeeweed removal subplots and a 1 m buffer around the subplots, all yankeeweed plants were removed by clipping stems and their resprouts in late 2017–early 2018. In 2017, before applying yankeeweed removal treatments, the presence/absence and percent cover of plant species in each subplot were recorded, and no significant differences in plant community composition between subplots were found (LaPierre, 2020). Likewise, no differences in plant community composition among larger, burn treatment plots were found (LaPierre, 2020).

Flower abundance and diversity

Flowers were sampled in each control and removal subplot per burn treatment plot in 2019 and 2020. Data collected in 2019 followed the first year of seasonal burn treatments (applied in 2018) and represented the year following fires in a typical biennial prescribed burn cycle, hereafter referred to as the year between prescribed burns or non-burn year (Fig. 1). Data collected in 2020 covered the second year following seasonal burn treatment initiation and represented a year in which fire was applied, hereafter referred to as the burn year (Fig. 1). Data were collected in all subplots once every 2 weeks, from early March to late November each year. Sampling every 2 weeks is frequent enough to capture differences in overall floral resource availability and flowering phenology across plant species (Platt, Evans & Davis, 1988). An observer recorded the flower density for plant species in bloom in each 10 m × 10 m subplot during each sampling event. For all non-composite flowering plant species, individual flowers were counted. Each flower head was counted as one unit or flower for composite flowering plant species (e.g., Asteraceae). For plant individuals with numerous (~1,000+) flowers, flower abundance was estimated by counting plant stems in bloom and multiplying by the average number of flowers per 10 stems per plant. Most plants were recorded to species (60/64), with a minority recorded to the genus level (4/64). Species richness of plants in bloom on each sampling date was determined as the number of species blooming per subplot.

Pollinator activity and diversity

Sampling for insect pollinators occurred from 8:00–13:00 when bees and other pollinators forage (Danforth et al., 2019). Observations only took place under partial to full sun, with no precipitation and temperatures 13 °C or greater. Due to the proximity of the treatment plots to one another, all within a 10-hectare area, we do not assume that treatment plots host independent pollinator communities given the large foraging range of many insect pollinators (Goulson & Stout, 2001; Pasquet et al., 2008). Thus, sampling was designed to measure differences in pollinator activity. To measure pollinator activity, pollinators were counted during 5-min timed observations per subplot on each sampling day in which flower density was also recorded. Two observers starting on opposite sides of the study area walked around each subplot in a slow, meandering manner for 5-mins per subplot while being careful not to disturb plants. In addition to standardized 5-min observation periods, observers also opportunistically recorded any insects seen on flowers during flower counts as described above. All insects were recorded to the following categories representing functional groups and/or broad taxonomic categories: honeybees (Apis mellifera), bumble bees (Bombus spp.), other bees (non-Apis, non-Bombus Apoidea, primarily solitary bees), wasps, flies (Diptera), butterflies and moths (Lepidoptera), and crawling insects (including ants and beetles). In addition to recording all pollinators to these broad categories, observers recorded the lowest identifiable taxonomic level for each insect. The order in which burn treatment plots were sampled changed across sampling days, and two observers working simultaneously ensured that all subplots within all burn treatment plots were sampled within a ~3-h period.

Flower visitors were not captured during the 5-min observation periods for the following reasons, which generally would have affected our ability to accurately assess the main response variable of interest, pollinator activity: (1) netting insects on low-lying plants is likely to disrupt other visitors in the area, (2) collecting insect visitors on the first observation of their foraging would limit latter observations of these same insects, (3) across such broad taxa of insects, some are much more readily captured than others, and (4) flower visitors were generally infrequent and thus lethal capturing was kept to a minimum. Instead, after standardized 5-min observation periods, in subplots where insect activity was observed, observers opportunistically collected a subset of insects from novel plant-pollinator interactions to identify them to the species level. As such, collecting did not have a standardized intensity but was instead done to gather descriptive information on plant-pollinator interactions in this ecosystem. All collected insects were pinned and identified using Discover Life, a reference collection of wild bees identified by entomologist Glenn Hall and taxonomist John Ascher, and an online key to bees of Florida (Pascarella & Hall, 2016).

Statistical analyses

We used generalized linear mixed effects models to examine the effects of timing of prescribed burn and yankeeweed removal treatment on flower density, flower richness, and pollinator activity (i.e., visitation rates). Models included the following fixed effects: seasonal burn treatment, yankeeweed removal treatment, day of year, interactions between seasonal burn treatment and day of year and between yankeeweed removal treatment and day of year, year (2019, 2020), interactions between seasonal burn treatment and year and between yankeeweed removal treatment and year, and time-since-burn (continuous, in days). In another model, we removed the density of yankeeweed from the overall flower density and used all fixed effects described above for other models, but also included a three-way interaction between seasonal burn treatment, yankeeweed removal treatment, and year to assess how the effect of yankeeweed treatment on flower density (excluding yankeeweed) varies across burn treatments and years. All models also included a random intercept of subplot (1–24) nested within plot (1–12) to account for correlated errors among subplots within plots repeatedly sampled within a year. We compared models with a Poisson distribution and a negative binomial distribution for each response variable and chose the best fit model (lowest AICc), which used a Poisson distribution for flower density and richness, and a negative binomial distribution for pollinator activity. Models were created with the “glmer” function using the lme4 package (Bates et al., 2014) in R. Statistical significance for each fixed effect was determined with Type II sum of squares ANOVA when interactions were not significant and Type III sum of squares when they were significant (see Lewsey, Gardiner & Gettinby, 2001; Langsrud, 2003).

We additionally examined year-long, cumulative richness of flowers using linear mixed models with the following fixed effects: seasonal burn treatment, yankeeweed removal treatment, year (2019, 2020), the interaction between seasonal burn treatment and year, and an interaction between yankeeweed removal treatment and year. Residuals of the year-long flower richness per subplot were normally distributed. Estimated marginal means and 95% confidence intervals for fire and yankeeweed treatments, and their respective interactions, were calculated using the effects package in R, and pairwise comparison of estimated marginal means were performed using the emmeans package (Lenth et al., 2021) in R. All residual diagnostics were computed using the DHARMa package (Hartig, 2018).

To examine the effects of treatments on flower community composition each year, we constructed a Bray-Curtis dissimilarity matrix using the year-long proportion of each flower species in each subplot. For this analysis, as the goal was to examine how yankeeweed removal influences the remaining community composition of floral resources, we removed yankeeweed from the community in each subplot. We calculated the proportions of the remaining total floral resources and ran a PERMANOVA on the Bray-Curtis dissimilarity matrix with the following explanatory variables: seasonal burn treatment, yankeeweed removal treatment, year (2019, 2020), the interaction between seasonal burn treatment and year, and an interaction between yankeeweed removal treatment and year. For significant effects, we followed the PERMANOVA with a SIMPER analysis to identify flower species that significantly contribute to differences between treatments. Analyses were done using the R package vegan (Oksanen et al., 2017).

Finally, we constructed plant-pollinator networks for each year using data combined across all plots and subplots using the R package igraph (Csardi & Nepusz, 2006). The network was created with weighted matrices containing interaction frequencies between each flowering plant species and pollinator taxa. For simplicity, we show all pollinator taxa but only plant taxa that received pollinator visits. Data informing these networks came from both timed observations and opportunistic pollinator collections across all plots and subplots.

Plant-pollinator interactions in this ecosystem

We recorded 64 and 69 flowering plant species in 2019 and 2020, respectively. Eupatorium compositifolium (yankeeweed), Dalea pinnata, Callicarpa americana, and Croton michauxii were the most abundant flowers in both years, although D. pinnata had fewer flowers in 2020. Insects visited a total of 39 and 46 flowering plant species across all subplots and sampling dates in 2019 and 2020, respectively (Appendix 1). In both years, insect visits were most frequent to the following plant species: C. americana, Chamaecrista fasciculata, Croptilon divaricatum, C. micauxii, D. pinnata, Rubus cunefolius, and Vaccinium stamineum. Insects visited more flowering plants in 2020, including the following additional species: E. compositifolium, Indigofera sp., Liatris sp., Solidago spp., Serenoa repens, and Stillingia sylvatica which received few to no pollinator visits in 2019. Dalea pinnata attracted honeybees, butterflies, and flies primarily, and C. americana attracted mostly bees and wasps. Additionally, E. compositifolium attracted primarily flies, and S. repens attracted mainly crawling visitors and flies. Other plants such as Liatris sp., C. fasciculata, C. michauxii, C. divaricatum, and Indigofera sp. attracted a variety of insects, including a diversity of bees (Appendix 1).

Effects of prescribed burn treatments on flowers and pollinators

Prescribed burn treatments had a significant effect on flower density, but this effect varied across sampling dates as indicated by a significant interaction (Table 1). Flower density generally increased throughout the season and peaked in October during wildflower bloom; differences in flower density across burn treatments also increased throughout the season and were most pronounced during peak bloom (Table 1). Furthermore, prescribed burns had varying effects on flower density across sampling years; in 2019 when plots were not burned (i.e., in between fire treatments), the summer-wet season fire treatment plots had the highest flower density, significantly higher than unburned control plots (Table 1; Fig. 2A), while in 2020 when plots were burned for a second time, the highest flower density was observed in the spring season fire treatment plots, significantly higher than unburned control plots (Table 1; Fig. 2A).

Flower richness did not vary significantly among burn treatments when analyzed per sampling date or as total richness per year (Table 1). Flower community composition also did not differ significantly among burn treatments (Table 1), and this was consistent across years (Table 1). However, flower communities did vary between years (Table 1). In 2019, D. pinnata, E. canadensis, and C. divaricutum were relatively more abundant, while in 2020, C. michauxii, S. repens, V. stamineum, C. americana, and Solidago spp. were relatively more abundant (Table 2).

The effect of prescribed burns on pollinator activity varied across years (Table 1; Fig. 2B) but was consistent across sampling dates within a year (Table 1). In 2019, pollinator activity was highest in the summer-wet and spring season fire treatment plots, both significantly higher than in unburned control plots (Table 1). In 2020 when fire treatments were again applied, pollinator activity was highest in the winter-dry and spring season burn treatments, both significantly higher than the unburned control plots (Table 1; Fig. 2B).

After accounting for the effects of prescribed burn treatments, time-since-fire had a positive effect on both flower density and flower richness as well as on pollinator activity (Table 1).

Effects of yankeeweed removal on flowers and pollinators

The effect of yankeeweed removal on flower density varied across sampling dates (Table 1) and years (Table 1; Fig. 3A). Overall, flower density was higher in control (no removal) plots in both years (Fig. 3A), and the blooming of yankeeweed largely drove this result. Yankeeweed removal did not affect flower richness per date (Table 1), total flower richness per year (Table 1), or flower community composition (Table 1). When we analyzed the effect of removal treatment on the flower density of all other plants (excluding yankeeweed), the effect of yankeeweed removal on flower density varied across burn treatments and years as indicated by the three-way interaction (Table 1; Fig. 4). In 2019, flower density was higher in yankeeweed removal plots compared to control plots, but only in the winter-dry season burn treatment (Fig. 4; Table 1). In 2020, flower density was less variable between yankeeweed removal and control plots (Fig. 4).

The effect of yankeeweed removal on pollinator activity varied across years (Table 1), with higher pollinator activity recorded in the yankeeweed removal plots in 2019 but comparable frequencies between treatments in 2020 (Fig. 3B).