Effects of Rhododendron removal and prescribed fire on bees and plants in the southern Appalachians

Abstract Rhododendron maximum is an evergreen shrub native to the Appalachian Mountains of North America that has expanded in recent decades due to past disturbances and land management. The purpose of this study was to explore how bees and plants were affected by the experimental removal of R. maximum followed by a prescribed fire in one watershed compared to a neighboring reference watershed. Bees and plants were sampled for three years in both watersheds. Comparisons were based on the rarefaction and extrapolation sampling curves of Hill numbers as well as multivariate methods to assess effects on community composition. Bee richness, Shannon's diversity, and Simpson's diversity did not differ between watersheds in the year after removal but were all significantly higher in the removal watershed in year two, following the prescribed fire. Bee Shannon's diversity and Simpson's diversity, but not richness, remained significantly higher in the removal watershed in the third year. Similar but weaker patterns were observed for plants. Comparisons of community composition found significant differences for bees in the second and third year and significant differences for plants in all three years. For both groups, significant indicator taxa were mostly associated with the removal watershed. Because bees appeared to respond more strongly to the prescribed fire than to the removal of R. maximum and these benefits weakened considerably one year after the fire, clearing R. maximum does not appear to dramatically improve pollinator habitat in the southern Appalachians. This conclusion is underscored by the fact that about one quarter of the bee species in our study area were observed visiting R. maximum flowers. The creation of open areas with wildflowers may be a better way to benefit bees in this region judging from the high diversity of bees captured in the small roadside clearings in this study.


| INTRODUC TI ON
Native bees and other pollinators within the historically forested regions of North America have experienced considerable changes since European colonization. Among these, the near-complete conversion of old growth forests to agriculture, development, or second-growth stands was particularly transformative (MacCleery, 1993). Studies exploring the relationship between surrounding forest cover and bee diversity have yielded mixed results. While some report that forest cover has a positive effect on the diversity of bees (Bennett & Isaacs, 2014;Taki et al., 2007;Watson et al., 2011), others suggest the opposite (Winfree et al., 2007). Such discrepancies may be driven by differences in habitat associations among bee taxa.
For example, Smith et al. (2021) found that about a third of all bees in the northeastern United States are associated with forests. Another third were open-habitat species while the remainder were classified as habitat generalists. Although previous studies from that region have shown a negative relationship between forest cover and total bee diversity (Winfree et al., 2007), Smith et al. (2021) reported a positive relationship between the amount of forest area and the diversity of forest-associated species. Moreover, Harrison et al. (2018) concluded that forest bees are largely replaced in agricultural and urban environments by species adapted to open habitats. Such findings suggest that forests play a key role in conserving bee diversity and underscore the importance of research aimed at optimizing forest conditions for these important and imperiled insects.
Today, forests in the United States cover about two thirds of the area that was forested in the year 1600 but differ from those earlier forests in many respects (MacCleery, 1993). The introduction of non-native species, fire suppression, fragmentation, and altered stand density have all changed forest composition and structure. It is apparent from a growing body of literature that native bees are sensitive to forest conditions (Hanula et al., 2016). One of the clearest patterns is that open stand conditions promote more diverse bee communities, in part through increasing floral resource availability near the forest floor (Hanula et al., 2016). As a consequence, more open early successional forests often support higher bee diversity than older forests with more closed canopies (Rivers & Betts, 2021;Taki et al., 2013). In mature forests, a number of management interventions have been shown to result in more open conditions and higher bee diversity. These include stand thinning, prescribed fire, and the mechanical removal of invasive shrubs (Hanula et al., 2016;. Invasive shrubs can include non-native species that dominate the midstory and largely exclude herbaceous plants from the forest floor as well as native species that have increased in abundance due to past disturbances and management history. One example of the latter concerns Rhododendron maximum L., an evergreen ericaceous shrub native to the Appalachian Mountains of the eastern United States. This shrub now dominates the forest midstory over large areas due to a combination of logging, the loss of American chestnut (Castanea dentata (Marsh.) Borkh.) and Eastern hemlock (Tsuga canadensis (L.) Carr.) from the overstory, and fire suppression Ford et al., 2012;Van Lear & Brose, 2002). Herein, we present the results from a study aimed at assessing the effects of R. maximum removal followed by prescribed fire on bee communities in North Carolina.
The Southern Appalachian region is considered a biodiversity hotspot, supporting North America's highest diversity of trees as well as many other plant and animal species (Stein et al., 2000). However, disturbances have altered forest structure and species composition of southern Appalachian watersheds (Adams, 2013;Bearup et al., 2014;Brantley et al., 2013Brantley et al., , 2015Ice et al., 2004). These include pre-and early settlement fires, logging in the early 20th century, drought, hurricanes, and insect and disease outbreaks. Most notable among these outbreaks was the introduced chestnut blight fungus (Endothia parasitica (Murr.) P.J. And. & H.W.) in the 1920s-1930s that led to widespread mortality of American chestnut. Consequently, overstory tree species composition changed from dominance by American chestnut in the 1930s, to drought-tolerant, xerophytic oaks (Quercus spp.) by mid-century, and finally to dominance by drought-intolerant, mesophytic species such as red maple (Acer rubrum L.) and tulip poplar (Liriodendron tulipifera L.) by the end of the century (Elliott & Swank, 2008;Elliott & Vose, 2011;Ford et al., 2012;Nelson, 1955). Subcanopy species composition changed along with the canopy changes, including the expansion of R. maximum.
Rhododendron maximum is an important species in these forests as it is highly shade tolerant, forms a dense subcanopy that shades the forest floor (Clinton, 2003), has little to no herbaceous or tree seedling (henceforth, understory) cover below its subcanopy (Beckage et al., 2000;Clinton et al., 1994), and it decreases N availability in the soil and litter layer to nonericaceous species (Wurzburger & Hendrick, 2007, 2009. As a result, tree height, biomass, and productivity are substantially reduced in areas with dense R. maximum cover (Bolstad et al., 2018). With regard to foraging bees and other pollinators, R. maximum flowers prolifically in the early summer months. Although the flower-visiting insect communities associated with other species of Rhododendron have been documented in previous studies from Europe and Asia (Sugiura, 2012;Tiedeken & Stout, 2015), no previous effort, to our knowledge, has been made to characterize the bee community visiting R. maximum in North America. This is a particularly interesting question given that the nectar produced by some Rhododendron species contains secondary compounds known to be toxic to some bee species (Tiedeken et al., 2016).
In this study, we used a paired-watershed study design to investigate the effects of R. maximum removal followed by prescribed fire on the diversity of bees and herbaceous plants. Based on the

T A X O N O M Y C L A S S I F I C A T I O N
Applied ecology; Biodiversity ecology; Community ecology; Conservation ecology; Ecosystem ecology; Entomology; Global change ecology; Invasion ecology; Restoration ecology literature summarized above, we hypothesized that the diversity of both groups would be significantly higher in the removal watershed and that the strongest effects would be detected after the prescribed fire. We also predicted that individual bee species would show positive associations with the removal watershed but not the reference watershed. We also sought to better understand the value of R. maximum to the local bee community by sampling bees directly from the flowers of this shrub.

| Study area and treatment
This study was conducted at the Coweeta Hydrologic Laboratory with increasing dominance of mesophytic species (e.g., red maple, tulip poplar) in the overstory and R. maximum in the subcanopy (Elliott & Swank, 2008;Elliott & Vose, 2011).
We utilized an ongoing study aimed at cutting R. maximum followed by prescribed fire as a way to increase water yield and understory plant diversity. The design follows the paired-watershed approach (Wilm, 1943) involving two neighboring watersheds with similar climate, vegetation, stand age, aspect, elevation, topography, and soils. The strength of this design is that the treatment effect can be isolated by accounting for year-to-year variation in climate using the reference watershed. The "removal watershed" is east-facing and mid-elevational (869-1141 m) with an area of 37.8 ha while the "reference watershed" is also east-facing and mid-elevational (920-1236 m) and 41.3 ha in area. Although neither watershed was uniformly covered by R. maximum at the beginning of the study, there was no significant difference in cover between the two watersheds initially based on observations at our 44 sampling plots (χ 2 = 0.76, p = .38). Another shrub, mountain laurel (Kalmia latifolia L.) was also present on the ridges in both watersheds and was not cut except in cases where it overlapped with R. maximum. As the Coweeta basin was last logged in the 1920s, before Coweeta Hydrologic Lab was established in 1934, both watersheds had been largely undisturbed for about a century or more prior to this study. with an aquatic label to avoid resprouting .
After the cut stems and slash had more than a year to dry, the same watershed was subjected to a prescribed fire on 23 March 2019 (Table 1). The reference watershed was not manipulated and served as the experimental control.

| Bee and plant sampling
We sampled bees and recorded herbaceous plant diversity in long-term permanent vegetation plots that were established in the 1930s. The plots within each watershed were arranged in two transects separated by ~400 m, and ~40 m separated our sampling locations within each transect. There were 23 and 21 plots in the reference and removal watersheds, respectively. To sample bees, we put out a set of three colored pan traps (white, yellow, and blue) in each plot, separated by 1 m and placed in a line parallel to For each sampling period, traps were collected after three days, pooled by plot, and returned to the lab for identification. We sampled in May, June, and August in 2018 and found bee numbers to drop quickly as the season progressed, a trend common for bees in temperate deciduous forests (Harrison et al., 2018). Therefore, we sampled monthly from April to June in 2019 and 2020 (Table 1).
All plants growing <0.5 m in height (regardless of flowering status) were recorded within a permanent 1 m 2 quadrat situated near the pan traps in late June to early July of each year (Table 1)

| Plot measurements
To better understand how local conditions influence bee and plant diversity, we collected data on canopy openness and burn severity at each plot. To measure canopy openness, we first used fisheye lenses to photograph overhead canopy cover. In 2018, these images were taken with a camera mounted on a tripod at ~1 m above the ground. Because the mountainous terrain made this approach cumbersome, we used a fisheye lens attachment on a cellphone (~1.7 m above the ground) to capture images in subsequent years.
Thus, canopy openness values from 2018 cannot be meaningfully compared to those taken in 2019 or 2020. We then used software (WinSCANOPY, regentsinstruments.com) to calculate canopy openness from each of these images. To obtain information on fire severity, we calculated the relative differenced normalized burn ratio (RdNBR, unitless) (Miller & Thode, 2007)

| Rhododendron maximum cover
Because neither watershed was uniformly covered in R. maximum, we recorded which of the plots in the removal watershed had been covered by the shrub pre-treatment based on the presence of cut R. maximum stems in the immediate vicinity of the bee sampling locations. For the reference watershed, we considered a plot to be partially covered by R. maximum if there was at least one stem >2.5 cm in diameter within a two-meter radius of any of our bowl locations or any R. maximum leaves or branches visible directly above any of the pan traps. We then assigned each plot to one of three R. maximum cover categories: removed (n = 10); absent (n = 19); and present (n = 15).

| Analysis
We first tested the abundance, richness, and Shannon's diversity of bees as well as plant richness and canopy openness for spatial autocorrelation. This was done by calculating Moran's I using the package ape (Paradis & Schliep, 2018) in R 3.6.1 (R Core Team, 2019). Because the variables were significantly autocorrelated (results not shown), and were thus not independent (an assumption of linear models), we used the iNEXT package (Hsieh et al., 2016) to compare bee and plant diversity estimates between the two watersheds based on the rarefaction and extrapolation sampling curves of Hill numbers. Hill numbers are a mathematically unified family of diversity indices (Chao et al., 2014) where the value of q determines how much weight is given to species based on their abundance. While q = 0 (richness) gives rare and abundant species equal weight, q = 1 (Shannon's diversity) gives more weight to abundant species and q = 2 (Simpson's diversity) counts only the dominant species (Hsieh et al., 2016

| Bees
In total, 110 species of bees were captured in this study (Appendices 3 and 4) of which a quarter were collected from R. maximum flowers ( Figure 2, Appendix 4) We found no differences in bee richness (q = 0), Shannon's diversity (q = 1), or Simpson's diversity (q = 2) in 2018, post-cutting but prior to the prescribed burn, based on overlapping confidence intervals at the base sample size (Figure 3). However, all three metrics were significantly higher in the removal watershed following the prescribed burn in 2019. These patterns were repeated in 2020 except for q = 0, which did not differ significantly between treatments ( Figure 3). Our separate comparison of R. maximum cover categories provides insights into the local effects of R. maximum cover on bees. There were no differences among the three categories in 2018 but plots from which R. maximum was absent or removed had significantly higher bee diversity at all levels of q than plots in which R. maximum was present in 2019 (Appendix 5). Except for no significant differences in species richness (q = 0), the results were similar in 2020. There was a tendency toward higher Simpson's diversity in the removed plots than in the absent plots, and this difference was significant in 2020. Bee richness increased with increasing plant richness in all three years (Table 2). In addition, bee richness increased with increasing canopy openness in 2019 and 2020 but not in 2018 (Table 2). Bee richness also increased with increasing burn ratio (RdNBR) in all three years in the removal watershed (Table 2). Our NMDS analyses failed to produce a useful ordination for bees in 2018 but yielded three-dimensional solutions for 2019 and 2020 (final stress = 16.74 and 16.97, respectively). Only the two axes explaining the most variation are presented in Figure 4, which for both years shows considerable overlap in bee communities between the two watersheds. According to MRPP, bee communities differed significantly between watersheds in the years following the prescribed fire, 2019 (T = −7.40, p < .001) and 2020 (T = −7.04, p < .001), but did not differ from one another after cutting, 2018 (T = −0.82, p = .19).
Indicator species analysis identified 19 bee species that were significantly associated with one of the two watersheds (Table 3). Eighteen were associated with the removal watershed beginning in 2019. Only one species, Bombus impatiens, was associated with the reference watershed and this was only the case in 2018.

| Plants
A total of 54 plant species were recorded in our herb plots over the course of the study (Appendix 6). Similar to the results for bees, we found no differences in plant richness (q = 0), Shannon's diversity (q = 1), or Simpson's diversity (q = 2) between watersheds in 2018 based on overlapping confidence intervals at the base sample size ( Figure 5). We only found significant differences between watersheds for q = 2 in 2019. In 2020, we still found no significant difference in plant richness between watersheds, but Shannon's diversity (q = 1) and Simpson's diversity (q = 2) were both significantly higher in the removal watershed ( Figure 5). Our separate comparison of R. maximum cover categories provides more insight into how plant diversity is affected by the presence R. maximum. There were no differences among the three categories in 2018 although it is noteworthy that the absent plots accumulated species more quickly than the removed plots (Appendix 7). In 2019, plant richness (q = 0) was significantly higher in absent plots than in present plots and Shannon's diversity (q = 1) and Simpson's diversity (q = 2) were both significantly higher in the absent and removed categories than in the present category. The results for 2020 were the same as those for 2019 except that Simpson's diversity only differed significantly F I G U R E 3 Rarefaction (solid lines) and extrapolation (dashed lines) of bee diversity in the reference and removal watersheds. Separate results are given for Hill numbers 0, 1, and 2. The results for species richness (q = 0) are shown in the left-most panels, whereas those for Shannon's diversity (q = 1) and Simpson's diversity (q = 2) are shown to the right.

F I G U R E 4
Ordinations from nonmetric multidimensional scaling showing differences in bee community composition among the sampling plots for 2019 and 2020 (note: a useful ordination was not found for 2018). Open and closed triangles represent removal and reference plots, respectively, while dots indicate individual bee species. The vector for bee richness is included in each ordination between the absent and present plots (Appendix 7). As stated above, as plant richness increased, bee richness increased; but, plant richness was unaffected by canopy openness. Regardless of year, plant richness was unaffected by burn ratio (RdNBR, data not shown).
According to MRPP, plant communities differed significantly be- identified 6 plant species that were significantly associated with one of the two watersheds (Table 3). Only one species, R. maximum, was associated with the reference watershed.

| DISCUSS ION
This study aimed to explore the effects of R. maximum removal followed by prescribed fire on bee and herb-layer plant communities.
Our results suggest that both taxa generally benefited from the treatment, although these effects did not become apparent until 2019 after the prescribed fire. While this delayed response may have been observed even without the prescribed burn, the conclusion that applying fire following shrub removal more strongly benefits bees than shrub removal alone is consistent with the findings of a previous study. Campbell et al. (2007) found flower-visiting insects to be significantly more abundant in plots that were both mechanically cleared of Kalmia and R. maximum in the shrub layer and subsequently burned than in plots that received just one or neither of these treatments. Many studies have reported significant increases in bee abundance and richness the year following a fire (Carbone et al., 2018;. Although this is often attributed to increased floral resource availability in recently burned areas, it could also simply be due to enhanced trap visibility against a charred background. The second explanation is more likely in this case considering our sampling took place soon after the burn in 2019, presumably too soon for differences in floral resource availability to come into play. Moreover, weaker differences in bee diversity between watersheds were observed in 2020, including no significant difference in species richness, after the plant community had had more time to recover from the burn. When analyzing data from the removal watershed separately, we found significant positive correlations between bee richness and fire severity (RdNBR) ( Table 2). Although this is consistent with previous studies (Galbraith et al., 2019), the fact that this was true even in 2018, the year before the burn, suggests that some correlate of fire severity, rather than fire severity itself, is likely responsible for these relationships. Because plant richness, canopy openness, and elevation were not significantly correlated with RdNBR, these metrics can be ruled out as factors driving these patterns. Although not measured in this study, it is possible that downed woody debris may have resulted in both higher local fire severity as well as higher bee richness. Dead wood represents one of the most important resources to forest insect communities (Ulyshen, 2018) and is known to provide critical nesting sites for several of the bee species captured in this study (e.g., Lasioglossum subviridatum (Cockerell), L.
The responses of plants to R. maximum removal and prescribed fire were similar to, but less dramatic than, those for bees. Although maximum cover classes show clearly that herbaceous plant diversity is higher in areas without R. maximum, especially soon after a fire.
Interestingly, by 2020 the rarefaction and extrapolation curves for plots from which R. maximum was absent or removed were nearly identical, suggesting a rapid recovery of herbaceous plant species richness following the removal of the shrub.
We found bee richness and plant richness to be significantly correlated in all three years. This is more likely due to both groups being more species rich in relatively open areas free from R. maximum than to bees responding directly to differences in plant richness. Many of the plant species recorded in our plots were either wind-pollinated species, too young to flower or were not observed flowering during the course of this study. Based on our general observations, floral resource availability is much lower within the closed forest than in clearings where we commonly observed bees actively visiting flowers. Although our sampling intensity was eleven-fold greater in the two watersheds combined than in the roadside clearings in 2019, the clearings yielded almost as many species (66 vs 67) and over a third as many individuals. These findings suggest that roadsides and other clearings may be particularly important foraging habitats for bees in the southern Appalachians and perhaps other regions dominated by mature closed-canopy forests.
We captured a diverse assemblage of bees visiting R. maximum flowers in this study, consisting of 27 species, and estimate that a quarter of all bee species in our study area utilize this floral resource.
Among these species is at least one species, Andrena cornelli Viereck, that specializes on R. maximum and closely related taxa. Although we did not sample bees visiting R. maximum blooms within the forest, it is clear from these results that R. maximum provides an important resource to bee communities in the southern Appalachians.

F I G U R E 5 Rarefaction (solid lines) and extrapolation (dashed lines) of plant diversity in the reference and removal watersheds.
Separate results are given for Hill numbers 0, 1, and 2. The results for species richness (q = 0) are shown in the left-most panels, whereas those for Shannon's diversity (q = 1) and Simpson's diversity (q = 2) are shown to the right. All curves include 95% confidence intervals, and comparisons are made at twice the smallest reference sample size (i.e., 42 sampling units for 2018 and 2020 and 38 units for 2019)

| CON CLUS IONS
Rhododendron maximum has benefited from the past century of disturbance in the southern Appalachians and is a more dominant member of the shrub layer than it was historically. The findings from this study suggest that both bees and herb-layer plant diversity will increase to some degree from local efforts to reduce R. maximum cover and reintroduce fire. The benefits of such measures to bees appear to be relatively mild and short term, however, especially compared to the dramatic increases in plant and bee diversity reported after the eradication of Ligustrum sinense Lohr, a non-native shrub . Most notably, we only detected a significant increase in bee species richness (q = 0) in 2019, the season immediately following the fire, whereas no sig-

CO N FLI C T O F I NTE R E S T
The authors declare there are no competing interests.   Nomada sp9* 4/3/2 0/5/3 10 27

A PPEN D I X 3 (Continued)
A PPEN D I X 5 Rarefaction (solid lines) and extrapolation (dashed lines) of bee diversity in plots from which R. maximum was absent, present, or removed.
Separate results are given for Hill numbers 0, 1, and 2. The results for species richness (q = 0) are shown in the left-most panels, whereas those for Shannon's diversity (q = 1) and Simpson's diversity (q = 2) are shown to the right. All curves include 95% confidence intervals, and comparisons are made at twice the smallest reference sample size (i.e., 20 sampling units)

A PPEN D I X 6
List of plant species recorded within the herb plots (within 0.5 m of the ground) in the reference and removal watersheds. The presence of a species is denoted with an x for each year Species Floral resource in herb layer?