Understory community shifts in response to repeated fire and fire surrogate treatments in the southern Appalachian Mountains, USA

Background: Decades of fire exclusion in the southern Appalachian Mountains, USA, has led to changing forest structure and species composition over time. Forest managers and scientists recognize this and are implementing silvicultural treatments to restore forest communities. In this study, conducted at the southern Appalachian Fire and Fire Surrogate Study site in Green River Game Land, North Carolina, USA, we assessed the effects of four fuelreduction methods (burned four times, B; mechanical treatment two times, M; mechanical treatment two times plus burned four times, MB; and control, C) on the changes in understory community from pre-treatment to posttreatment years (2001 to 2016). We used non-metric multidimensional scaling (NMDS) to determine overall understory community heterogeneity, agglomerative hierarchical cluster analyses (AHCA) to determine finer-scale changes in understory community structure, and indicator species analyses (ISA) to identify the species that were associated with the different fuel reduction treatments over time. Results: The NMDS ordination showed little separation between treatment polygons. The AHCA resulted in two main categories of understory species responses based on how treatment plots clustered together: (1) species apparently unaffected by the treatments (i.e., no treatment pattern present within cluster); and (2) species that responded to B, M, or MB treatments (i.e., pattern of treatment plots present within cluster). Nearly half (49.2%) of tree-species plots clustered based on treatments; 60% of shrub-species plots clustered based on treatments; and 64% of herbaceous-species plots clustered based on treatments. Many plots clustered similarly in response to firerelated treatments (B and MB). The ISA identified 11 total tree species: three in B, one in M, and seven in MB; six total shrub species: two in M, and four in MB, and 17 total herbaceous species or genera: one in C, and 16 in MB. Conclusion: Fire and fire surrogate treatments did not dramatically shift understory composition after 15 years. However, certain ruderal and early seral species responded positively to MB, which was the most intensive treatment. Modest understory community changes were also observed in B, suggestive of early signs of shifting composition toward a more open forest community after four burns.


Background
Natural and anthropogenic disturbance regimes influence forest ecosystems by altering resource distribution, which changes vegetation structure and composition across a range of spatial and temporal scales (Pickett et al. 1999;Certini 2005). Variations in disturbance type, intensity, and frequency often create mosaics in a landscape, increasing community-level heterogeneity (Willig and Walker 1999;Greenberg et al. 2016). Repeated disturbance has shaped many complex ecosystems over time, often resulting in the development of disturbance-dependent forest communities (Horn 1974). Historical records suggest that the consistent presence of fire in a forested landscape favors the establishment, growth, and dominance of certain plant species (Greenberg et al. 2016;Lafon et al. 2017). Thus, fire acts as a selective pressure that has driven the evolution of adaptations and physiological strategies that make them more resilient to-if not dependent upon-burning (Pyne 1997). Over time, these traits have allowed certain species to remain dominant in many fire-maintained forest communities. Since fire suppression began in the early 1900s, many forest ecosystems throughout the US have experienced increased fuel loading, increased wildfire risk, changes in vegetation composition, and landscapelevel homogenization (Huddle and Pallardy 1999;Taylor 2007). In the United States, one of the regions most heavily affected by fire suppression is the southern Appalachian Mountains (Garren 1943;Greenberg et al. 2016).
Characterized by extreme topography, varying soil types, and large precipitation gradients, the southern Appalachian Mountain region supports high levels of biodiversity and heterogeneity (Van Lear and Waldrop 1989;Waldrop et al. 2014). However, these complexitiesin concert with extensive exurban development-also make the southern Appalachian Mountains one of the most difficult regions to manage with fire (Stanturf et al. 2002). Due to fire suppression, forests in the southern Appalachian Mountains have experienced increases in midstory and overstory density, excess fuel loading, oak (Quercus L.) and yellow pine (Pinus L. subgenus Diploxylon) regeneration failure, forest homogenization, and encroachment of mesic hardwood species (Brose and Van Lear 1998;Elliott et al. 1999). A prevalent structural shift from open pine−hardwood forests and woodlands to predominantly closed-canopy mesic forests has been documented in this region (Waldrop et al. 2008). These open forests and woodlands were typically composed of firedependent species, such as oaks, yellow pines, and diverse herbaceous flora (sedges, grasses, and forbs; Ayers and Ashe 1905). Conversely, contemporary mesic hardwood forests support the growth of fire-sensitive species like red maple (Acer rubrum L.) and yellow-poplar (Liriodendron tulipifera L.) (Nowacki and Abrams 2008). Additionally, there have been increases in flammable midstory shrubs like mountain laurel (Kalmia latifolia L.) and rhododendron (Rhododendron L. spp.), which outcompete firedependent species and create potential for higher intensity fires (Monk et al. 1985;Brose et al. 2001).
Despite a majority of vascular plant diversity being found in the understory layer, this forest component is comparatively understudied in southern Appalachian Mountain forests (Gilliam and Roberts 2003). Understory plant communities in this region can be composed of a wide array of tree, shrub, and herbaceous species regeneration (Waldrop et al. 2016;Oakman et al. 2019). This understory community is largely shaped by topography-related microsite characteristics that include abiotic factors such as microclimate, soil moisture, and soil fertility (Hutchinson et al. 1999). Repeated fires of varied intensity and severity often contribute to fluctuations in these microsite characteristics by increasing light availability from canopy openings, which alters understory community composition (Small and McCarthy 2002;Gilliam and Roberts 2003). However, many biotic factors also contribute to the formation of understory community composition (Azeria et al. 2011). For example, species-specific characteristics, like shade intolerance or prolific resprouting abilities, can influence post-fire succession (Van Lear et al. 2000;Hutchinson et al. 2005).
Many herbaceous species are classified as early seral and ephemeral, and thus respond positively-albeit temporarily-to canopy recently opened by disturbance (Roberts 2004;Hutchinson et al. 2005). These species often grow well in exposed site conditions like those created by canopy disturbance (Pavlovic et al. 2011). Understory shrubs often respond positively to cutting or dormant-season burning by sprouting after topkill (Chapman 1950;Moser et al. 1996). Competition from these shade-tolerant shrubs often impedes the regeneration of oaks, yellow pines, and herbaceous species (Monk et al. 1985). Understory hardwoods and softwoods often respond positively under open conditions in the mid-and overstory created by disturbance (Kuddes-Fischer and Arthur 2002). For example, hardwoods like oaks and softwoods like yellow pines are more tolerant of fire and respond more positively to higher intensity fire treatments than mesic hardwoods like red maple and mesic softwoods such as eastern white pine (Pinus strobus L.), particularly in drier sites (Elliott and Vose 2005;Holzmueller et al. 2014). Mesic hardwood or softwood presence reflects an increase in canopy density, as upland oaks typically regenerate well on more xeric or open sites; most mesic hardwoods are shade-tolerant or generalist species and regenerate well in a variety of site types (Huddle and Pallardy 1999;Iverson et al. 2008). Characteristics like shade and moisture tolerance influence species-level interactions, often causing shifts in species composition (Blankenship and Arthur 2006). Therefore, identifying these species traits may inform our understanding of how vegetation will respond to different fuel-reduction treatments (Keyser et al. 2008;Azeria et al. 2011).
In the southern Appalachian Mountains, prescribed fire is used for objectives such as reducing wildfire risk, creating early seral wildlife habitat, and restoring fire-adapted plant communities (Monk et al. 1985;Christensen 1993;Waldrop et al. 2016). Although the most prevalent technique for achieving these objectives in the southern Appalachian Mountains has been dormant-season burning (January through March), scientists and managers in the region have also expressed interest in using other fuel reduction methods such as mechanical treatments like mastication or chainsaw felling of ladder fuels to reverse the effects of fire suppression (Brose and Van Lear 1998;Schwilk et al. 2009). Alone or in combination with prescribed fire, these treatments have shown promise for reducing fuels and restoring a more open forest structure, but some knowledge gaps remain (Waldrop et al. 2016). These include effects on oak and pine regeneration, the control of competitive fire-sensitive species, and the effects on herbaceous vegetation (Stephens et al. 2012;Waldrop et al. 2016).
The objective of this study was to investigate changes in understory plant communities in the southern Appalachian Mountains in response to three fuel-reduction treatments applied over the course of 15 years: burn only (B), mechanical (M), and a combination of mechanical and burning (MB). Unburned controls (C) were also included. We used ordination and plot-level clustering techniques to characterize how understory vegetation responds to the different treatments. Additionally, we used an indicator species analysis to identify species that were most strongly associated with the treatments. We hypothesized that, after 15 years, there would be a distinctive vegetation "signature" for each treatment, reflected in shifts in their ordination polygons, and vegetation similarities causing plots to cluster together by treatment. Further, we hypothesized that each treatment would contain certain indicator species-species common in that treatment but rare in others-and for the treatments that involved fire (MB and B), ruderal, early seral, and fire-adapted species would be prevalent among them.

Study area
The study area was located in Polk County, North Carolina, USA, on the Green River Game Land, which is managed for wildlife, public recreation, timber, and other resources by the North Carolina Wildlife Resources Commission (NCWRC; Fig. 1). The Green River Game Land covers 5841 ha and is classified as a mountainous region, where elevations range from about 300 to 800 m. Forests in the region are typical of the Southern and South-Central Oak-Pine Forest and Woodland Ecosystem macrogroup (NatureServe 2020). When the study was initiated in 2001, forests in the study area were about 80 to 120 years old and consisted of mixed xeric or mesic oak and pine species depending on the topographic position. Shortleaf pine (Pinus echinata Mill.), pitch pine (Pinus rigida Mill.) and Virginia pine (Pinus virginiana Mill.) can be found on dry ridge tops while eastern white pine (Pinus strobus) was found in moist coves. Ericaceous shrubs like mountain laurel (Kalmia latifolia) and great rhododendron (Rhododendron maximum L.), composed a dense midstory layer throughout the study area, with the former more prevalent in xeric sites and the latter more prevalent in mesic sites. Most of the soils are of the Evard series (fine-loamy, oxidic, mesic, Typic Hapludults) in areas that can be described as moderately deep, welldrained, mountain uplands. Bedrock is composed of igneous and high-grade metamorphic rocks, including mica gneiss, hornblende gneiss, and amphibolite (Web Soil Survey 2020).

Study design
This study utilized a randomized complete block design, with four treatment units in each of three replicate blocks for a total of 12 treatment sites. Each treatment site covered an average of 12 ha, which included a 4 ha buffer zone. Both the buffer zones and the treatment sites received the same experimental treatment, but no sampling plots were located within the buffer zones. Within the replicate blocks, each of the four treatment units were randomly assigned to one of the treatments: control (C), prescribed burning only (B), mechanical fuel reduction (M), and prescribed burning plus mechanical fuel reduction (MB). Treatment units were designed to include all prevailing combinations of elevation, aspect, and slope. However, these conditions varied within treatment units and were not separated for analysis (Fig. 1).
The B treatment was repeated four times by 2016, having been applied in February or March of 2003March of , 2006March of , 2012March of , and 2015. All fires were burned with a spot fire technique; the first was done by helicopter ignition and the others were done by hand ignition. Fire intensity was generally low, with flame lengths typically <1 m. Localized areas of higher fire intensity (flame lengths >2 m) were occasionally observed in MB, presumably due to higher fuel loading created by that treatment.  (2001), and in the growing seasons following each treatment (2003, 2005, 2006, 2008, 2011, 2012, 2013, 2014, 2015, and 2016). Visual estimates indicated nearly 100% burn coverage in both B and MB treatments during each fire.
A 50 × 50 m grid was established in the treatment areas, with grid points permanently marked and georeferenced. Ten 0.1 ha sample plots were established at randomly selected grid points within each treatment area (Fig. 2). The sample plots were 50 × 20 m and divided into ten subplots, each about 10 m 2 . Within each subplot, two 1 m 2 quadrats were established in the northwest and southeast corners to measure understory vegetation (<1.4 m tall) using Modified Whittaker plots (Oakman et al. 2019). Composition and abundance data were collected in each 1 m 2 quadrat; this included recording the cover values of all species and additionally recording the stem counts for tree species. Cover values used in this sampling method were visually estimated and recorded as classification values: 1 (<1%), 2 (1 to 10%), 3 (>10 to 25%), 4 (>25 to 50%), 5 (>50 to 75%), and 6+ (>75%). To generate workable cover class values for analysis, we used the class midpoint of the percent ranges for each cover class; for example, 5.5% would be used for the cover class 2, etc.
We acknowledge that the asynchronous nature of the experimental treatments potentially confounds our This asynchrony is an unavoidable reality of applying landscape-level treatments; it would have been logistically and financially challenging to conduct more than two mechanical treatments during the study period.

Analysis
All species in the understory (<1.4 m tall) were classified into three general life form categories for analysis: trees, shrubs, and herbaceous species. Woody vines were excluded to reduce bias in the dataset, as collection methods were inconsistent in 2001 and 2016. Cover for trees, shrubs, and herbaceous species were derived from cover classification values recorded in the posttreatment year, 2016. The midpoint percentage represented in each cover class was then averaged across each plot for all treatments (n = 120; see Oakman et al. 2019 for details). Stems per hectare (stem density) of each tree species was also calculated for each plot.
To assess understory community response to different long-term repeated fuel-reduction treatments, we first visualized the data in a non-metric multidimensional scaling (NMDS) ordination (Kruskal 1964a, b;McCune et al. 2000). Cover data from the 2016 survey period was used to construct a Bray-Curtis dissimilarity matrix for all understory species (Faith et al. 1987;Oksanen et al. 2015). Each life form category was analyzed using the metaMDS function from the vegan package in R software and overlaid in a single figure (version 3.2.2, R Core Team 2016). Species that occurred in less than 5% of all plots were excluded from the analysis. Stability for each life form category was assessed using the scree.plot function with 20 randomized runs and three final axes iterations. A two-dimensional solution was used in the final figure, as it was found to be sufficient for explaining similarities among species along ecological gradients. The final ordination contained three components: maximum convex polygons, species points, and axes. The convex polygons were associated with each treatment, representing variation in understory plant community responses and heterogeneity among plots in each treatment. The species points were positioned within the ordination to reflect treatments with which they are most closely associated. Axes can represent complex environmental gradients that describe separation among species points and treatment polygons.
We then conducted agglomerative hierarchical cluster analyses (AHCA), which form clusters based on shared species within our plots (McCune and Grace 2002). The cluster analyses were performed on a finer scale to show understory community shifts that may have otherwise been hidden within the NMDS ordination, allowing us to examine defined groups of species without the influence of environmental gradients. AHCAs were performed with the agnes function in the cluster package (Maechler et al. 2015). Each of the three life form categories were analyzed separately at the plot level (n trees = 120, n shrubs = 120, n herbaceous species = 117) to compare between-plot similarity of plant composition for different life form categories (Oksanen et al. 2015). Three plots were omitted within the herbaceous species dataset due to the absence of understory vegetation in those plots. Stem density data were used in the tree species AHCA, and cover data were used in the shrub and herbaceous species AHCA. Cover was normalized, making the marginal sum of squares equal to zero for a better fitting distribution. With an AHCA technique, each plot is considered an individual cluster and therefore plots with similar species are grouped into larger clusters resulting in a single dendrogram (McCune and Grace 2002). We used Bray-Curtis as the distance metric and flexible beta linkage (β = −0.25) as the fusion strategy to determine the appropriate number of clusters for each life form category based on local group structure (Oksanen et al. 2015). The cluster number was determined from fusion height, a visualization method that shows natural breaks in the data, indicating the highest number of plot similarities. To better explain plant community responses to the four treatments, clusters were then grouped into two descriptive categories based on the proportions of treatment plots within each cluster: (1) species apparently unaffected by the treatments, and (2) species that responded to treatments (McCune et al. 2000;Gonzalez-Tagle et al. 2008).
Indicator species analyses (ISA) were conducted in the indicspecies package in R for each life form category to identify indicators of species composition shifts in response to our treatments (Dufrêne and Legendre 1997;De Cáceres and Legendre 2009). Indicator species are those that have had specificity and fidelity to a given treatment and, thus, are the most indicative members of that treatment (Costanza et al. 2017). For shrub and herbaceous species, cover data were used to calculate the highest indicator value (IVmax), P-value, specificity, and fidelity for each species using the multipatt function (duleg = TRUE) and strassoc function in R (De Cáceres and Legendre 2009). A species' indicator value (IV; 0 to 1) is the square root of the product of that species' A and B; A indicates the probability that the given species is in a given cluster when it is found, and B indicates the probability of finding the given species in a given cluster (Shearman et al. 2017). The P-values (α = 0.05) represent the probability of obtaining an indicator value by chance that is equal to 1 (Kane et al. 2010

Non-metric multidimensional scaling (NMDS)
The NMDS ordination was resolved by two axes, with a stress value of 0.20, which is considered acceptable based on Clarke (1993, Fig. 3). Treatment polygons overlapped considerably within the ordination space, with relatively low distances between species. Overall, the resulting NMDS showed little separation between treatment polygons and little variability in species spread.

Agglomerative hierarchical cluster analysis (AHCA)
Considering the NMDS results, we used an AHCA to further break down vegetation life forms into discrete clusters at the plot level (P). Clustering plots with similar composition helped to describe small-scale patterns in plant community composition that may have otherwise been overlooked within the NMDS ordination (Shearman et al. 2017). These cluster-level responses were generalized into two broad categories based on species composition within the plots: (1) species apparently unaffected by the treatments (i.e., no apparent pattern in relationship to treatments; similar to C), and (2) species that responded to treatments (i.e., apparent pattern associated with treatments; different from C).
Among the clusters that responded differently to C, one or more distinct responses can be described by treatment-plot proportions (e.g., if a cluster contains B or MB plots). Clusters 2 and 5 in the tree dendrogram were associated with C (51% of total plots [TP Total ]), while clusters 1, 3, and 4 responded differently to B, M, and MB treatments (49.2% TP Total ) (Fig. 4). Plots in cluster 1 (T C1 ) appeared to respond similarly to the effects of B and MB (T C1 = 11 plots in B, six in C, two in M, and 15 in MB), plots in cluster 3 (T C3 ) responded similarly to the effects of B (six plots in B, one in C, one in M, and 0 in MB), and plots in cluster 4 (T C4 ) responded distinctly to M effects (zero plots in B, seven in C, ten in M, and zero in MB).
In the shrub dendrogram, cluster 4 responded similarly to C (16% of total plots [SP Total ]), while clusters 1, 2, and 3 responded differently to B, M, and MB treatments (60% SP Total ) (Fig. 5). Plots in clusters 1 (S C1 ) and 2 (S C2 ) responded similarly to the effects of B and MB (S C1 = six plots in B, three in C, one in M, and seven in MB; S C2 = ten plots in B, six in C, four in M, and nine in MB), and plots in cluster 3 (S C3 ) responded distinctly  A Bray-Curtis approach was used as the distance metric and a flexible beta linkage was used as the fusion strategy to determine the appropriate number of clusters. The cluster number was determined from fusion height, a visualization method that shows natural breaks in the data, indicating the highest number of plot similarities. Each line on the dendrogram is denoted by orange (burned four times; B), black (control; C), purple (mechanical treatment two times; M), and blue (mechanical treatment two times plus burned four times; MB) dots that indicate which treatment was applied to that plot. Cluster 1 had 11 plots in B, six in C, two in M, and 15 in MB, categorizing it as having a different response from C (category 2; Cat. 2). Cluster 2 had 13 plots in B, 15 in C, 16 in M, and 14 in MB, categorizing it as having a similar response to C (category 1; Cat. 1). Cluster 3 had six plots in B, one in C, one in M, and 0 in MB, falling under category 2. Cluster 4 had zero plots in B, seven in C, ten in M, and zero in MB, falling under category 2. Cluster 5 had zero plots in B, one in C, one in M, and one in MB, falling under category 1. The horizontal axis at the bottom of the dendrogram represents the distance or dissimilarity between clusters A Bray-Curtis approach was used as the distance metric and a flexible beta linkage was used as the fusion strategy to determine the appropriate number of clusters. The cluster number was determined from fusion height, a visualization method that shows natural breaks in the data, indicating the highest number of plot similarities. Each line on the dendrogram is denoted by orange (burned four times; B), black (control; C), purple (mechanical treatment two times; M), and blue (mechanical treatment two times plus burned four times; MB) dots that indicate which treatment was applied to that plot. Cluster 1 had six plots in B, three in C, one in M, and seven in MB, categorizing it as having a different response from C (category 2; Cat. 2). Cluster 2 had ten plots in B, six in C, four in M, and nine in MB, falling under category 2. Cluster 3 had eight plots in B, 15 in C, 23 in M, and nine in MB, falling under category 2. Cluster 4 had six plots in B, six in C, two in M, and five in MB, categorizing it as having a similar response to C (category 1; Cat. 1). The horizontal axis at the bottom of the dendrogram represents the distance or dissimilarity between clusters Fig. 6 A dendrogram of clustered plots (n = 117) derived from the agglomerative hierarchical cluster analysis (AHCA) for herbaceous vegetation shows five distinct clusters (labeled 1 through 5) in response to treatments (control = C, burned four times = B, mechanical removal two times of stems <10 cm dbh = M, and combination of M and B = MB). Data from 2016 were used in this study (2001 to 2016) of plant community change in response to fuel treatments in Green River Game Land, North Carolina, USA. A Bray-Curtis approach was used as the distance metric and a flexible beta linkage was used as the fusion strategy to determine the appropriate number of clusters. The cluster number was determined from fusion height, a visualization method that shows natural breaks in the data, indicating the highest number of plot similarities. Each line on the dendrogram is denoted by orange (burned four times; B), black (control; C), purple (mechanical treatment two times; M), and blue (mechanical treatment two times plus burned four times; MB) dots that indicate which treatment was applied to that plot. Cluster 1 had nine plots in B, 14 in C, 12 in M, and zero in MB, categorizing it as having a different response from C (category 2; Cat. 2). Cluster 2 had eight plots in B, zero in C, zero in M, and 20 in MB, falling under category 2. Cluster 3 had one plot in B, six in C, four in M, and one in MB, falling under category 2. Cluster 4 had five plots in B, five in C, five in M, and two in MB, categorizing it as having a similar response to C (category 1; Cat. 1). Cluster 5 had seven plots in B, three in C, eight in M, and seven in MB, falling under category 1. The horizontal axis at the bottom of the dendrogram represents the distance or dissimilarity between clusters Table 1 Site specificity, fidelity, relative indicator value (IV rel ; specificity × fidelity), abundance (stems ha −1 or percent cover), and the maximum indicator value (IV max ) results from the indicator species analyses (ISA) on understory tree, shrub, and herbaceous species in our study (2001 to 2016) of plant community change in response to fuel treatments in Green River Game Land, North Carolina, USA. These results, based on 2016 data from the study, indicate species compositions in response to repeated treatments (control = C, burned four times = B, mechanical removal two times of stems <10 cm dbh = M, and combination of M and B = MB). Asterisks represent the degree of significance (*** = P < 0.001, ** = P < 0.01, * = P < 0.05) In the herbaceous species dendrogram, clusters 4 and 5 responded similarly to C (36% of total plots [HP Total ]), while clusters 1, 2, and 3 responded distinctly to MB, B, and M treatments (64% HP Total ) (Fig. 6) (Table 1).

Non-metric multidimensional scaling (NMDS)
Contrary to our hypothesis, the overlapping treatment polygons suggest similarity in post-treatment species composition because treatments share many species. Most of the tree, shrub, and herbaceous species showed little separation in ordination space, suggesting somewhat differential but mostly shared responses to treatments. Additionally, the relatively short distances between species made the ecological trends represented by each axis difficult to interpret. Little separation between treatment polygons and little variability in species spread was also observed when all species were included in the ordination. This suggests that there may not have been a strong species response to the treatments and that further analyses may be necessary to examine the changes that may have occurred.

Agglomerative hierarchical cluster analysis (AHCA)
While the NMDS showed only modest differences in understory community composition between treatments, the results of the AHCA supported our hypothesis, showing evidence of treatment effects within life form categories. The largest portion of clusters was associated with category 1 (those with no apparent pattern in relationship to treatments; similar to C), suggestive of the predominance of generalist species that are largely unaffected by treatment type. However, some clusters responded similarly to fire-related disturbance (B or MB), suggestive of modest compositional shifts to a more ruderal or early seral plant community in those treatments. This may predominantly be due to the overstory and midstory canopy openings that were mostly created by MB treatments (due to localized areas of higher-intensity fire in MB), as reported in Waldrop et al. (2016). However, identifying individual species that are driving these responses would give more indication of compositional and abiotic changes in response to repeated treatments (Keyser et al. 2008;Azeria et al. 2011). Overall, many clusters showed similarities in vegetation response among treatments; however, the few clusters that showed divergence from the C treatments suggest only subtle and perhaps localized treatment effects on understory vegetation.

Indicator species analysis (ISA)
Within the tree group, the indicator species in B (Acer rubrum, Amelanchier arborea, and Liriodendron tulipifera) indicated that B sites were largely differentiated by mesic species that grow well under conditions created by lowintensity prescribed fire. Four repeated dormant-season burns did not appear to be sufficient for meeting the management objective of creating understory conditions that favor oak and yellow pine recruitment (Kuddes-Fischer and Arthur 2002; Dolan and Parker 2004). However, it is possible that fires of greater intensity or in a different season could produce a different result. The indicator species in M (Pinus strobus) suggests that these sites were differentiated by white pine, a mesic species (Phillips et al. 2007). This also suggests that long-term M treatments may not create conditions that are favorable for fire-tolerant species as they need a more open canopy and drier microsite conditions to grow optimally (Vose et al. 1993). The indicator species in MB (Sassafras albidum, Diospyros virginiana, Nyssa sylvatica, Oxydendrum arboreum, Quercus coccinea, Q. montana, and Robinia pseudoacacia) indicate that these sites are differentiated by more xeric species, many of which are light-responsive and grow well under more open conditions (Clinton and Vose 2000). More specifically, Quercus montana, Q. coccinea, and Sassafras albidum and Robinia pseudoacacia grow best in open, dry conditions, suggestive of some level of mesophication reversal in MB (Boring and Swank 1984;Dey and Hartman 2005), perhaps caused by greater fire intensity in that treatment.
Within the shrub group, the indicator species in M (Kalmia latifolia and Rhododendron maximum) indicate that these sites are differentiated by ericaceous shrubs that grow well in shade, prefer mesic conditions (R. maximum), and resprout prolifically when cut (Vose et al. 1993). This suggests that M treatments may not have reduced ericaceous shrub competitors, which is a priority management objective in this region (Waldrop et al. 2016). The indicator species in MB (Ceanothus americanus, Hypericum hypericoides, Lyonia ligustrina, and Rhus glabra) suggests that these sites are differentiated by more light-responsive and opportunistic species that grow well in open, xeric sites (Hutchinson et al. 2005;Keyser et al. 2008).
Within the herbaceous group, the indicator species in C (Arundinaria appalachiana) was widely prevalent, whereas it was sparse in other treatments. This suggests that Arundinaria appalachiana may be sensitive to disturbance, or that it is a poor competitor in a postdisturbance environment. The indicator species in MB suggest a high level of understory diversity, with five graminoids (Dichanthelium spp., Piptochaetium avenaceum, Schizachyrium scoparium, Carex sp., and Scleria spp.), and 11 forb species (Cassia sp., Conyza canadensis, Coreopsis major, Desmodium nudiflorum, Erechtites hieraciifolius, Helianthus divaricatus, Houstonia purpurea, Lespedeza bicolor, Potentilla canadensis, Rubus argutus, and Solidago sp.), three of which are nitrogen-fixing (Cassia sp., Desmodium nudiflorum, and Lespedeza bicolor). This suggests that MB treatments are facilitating the establishment of a different set of species that are largely unique to this treatment (Burton et al. 2011). The indicator species in MB are also suggestive of a shift towards annual or early seral species, which is a desirable outcome for fire managers (Waldrop et al. 2016). The forbs, mainly in the Asteraceae family, often respond well to larger disturbances, indicating larger openings in the canopy and midstory (Hutchinson et al. 2005). Many graminoids, such as Schizachyrium scoparium, often grow well in disturbed sites with high light availability (Peterson et al. 2007).

Conclusions
Four fires and two mechanical treatments over the course of 15 years resulted in only modest changes in understory vegetation composition, as observed from our NMDS results. Nonetheless, we observed the greatest degree of change, including an increase in early seral, fire-adapted, or fire-dependent understory species in the most intensive treatment (MB). This treatment likely resulted in the greatest increases in understory light availability, as well as reductions in litter and duff necessary for the establishment of these species. Further research, which should include the continued frequent application of prescribed fire, should be conducted on the longer-term effects of B to determine if the effects of B will eventually approach those of MB. Additionally, the differences in understory responses observed between M and B treatments suggests that M is only somewhat of a surrogate for B. The results of this study will prove valuable for managers in the southern Appalachian Mountain region who are considering using fire and fire surrogate treatments to manipulate vegetation composition.
Additional file 1. Species codes and scientific names for all taxa shown in the non-metric multidimensional scaling figure (Fig. 3). Data from 2016 were used in this study (2001 to 2016) of plant community change in response to fuel treatments in Green River Game Land, North Carolina, USA.

Acknowledgements
In addition to the funding sources listed above, we would like to acknowledge the current and past students, researchers, and agency employees that established the study plots, measured vegetation, and assisted with data management and analysis. We are grateful to each, but special thanks go to M. Akridge, G. Chapman, E. Gambrell, S. McKinney, H. Mohr, R. Phillips, D. Simon, and M. Smith. Feedback from two anonymous reviewers greatly improved the original manuscript. The manuscript is technical contribution number 6937 of the Clemson University Experiment Station. This material is based upon work supported by USDA National Institute of Food and Agriculture, under project number SC-1700514.
Authors' contributions ECO conducted the data analyses and wrote the first draft of the manuscript. TAW conceived of the study, secured initial and continued funding, and provided pre-fire data from 2001. DLH, TAW, and KB assisted with data analysis and edited the manuscript. All authors read and approved the final manuscript.

Funding
This study was made possible by generous funding and support from the National Fire Plan (NFP), Joint Fire Science Program (JFSP), US Forest Service, and the North Carolina Wildlife Resources Commission. The NFP and JFSP provided funding for initial plot establishment and data collection, and the US Forest Service provided funding for data collection in 2016 (Project #15-JV-11330136-044). The North Carolina Wildlife Resources Commission conducted the prescribed burns and assisted with general project management. This paper is based upon work supported by the USDA Hatch program, through the Clemson Experiment Station, under project number SC-1700514.
Availability of data and materials Please contact author for data requests.