Drivers of plant functional group richness and beta diversity in fire‐dependent pine savannas

To assess the drivers of plant functional group richness and beta diversity in fire‐maintained North American Coastal Plain (NACP) savannas.


| INTRODUC TI ON
The North American Coastal Plain (NACP), a global biodiversity hotspot due to its high vascular plant diversity and endemism, is a fire-prone region in which conservation depends upon the widespread application of prescribed fire to maintain native savanna and grassland biodiversity (Noss et al., 2015). At last count, over three million hectares of NACP savanna were in active prescribed fire rotation (Mitchell & Duncan, 2009), representing a massive ecosystem management achievement following decades of anthropogenic fire suppression (Van Lear, Carroll, Kapeluck, & Johnson, 2005). The reintroduction of fire has yielded a regional landscape dotted with fire-maintained savanna fragments. Many studies have provided valuable information about the factors affecting within-preserve plant species diversity and composition (i.e., fire frequency, land use history, hydrologic and soil gradients, and overstory-understory relationships) (Glitzenstein, Streng, Masters, Robertson, & Hermann, 2012;Kirkman, Mitchell, Helton, & Drew, 2001;Platt, Carr, Reilly, & Fahr, 2006;Veldman et al., 2014), but data on plant diversity patterns at the regional metacommunity level are lacking. With so many savanna preserves in the NACP now returned to a fire-maintained grass-tree structure, expanding to metacommunity-level biodiversity questions is timely.
Regional approaches to biodiversity conservation are increasingly recommended for stemming global biodiversity loss in the face of accelerating anthropogenic change (Socolar, Gilroy, Kunin, & Edwards, 2016). Fire-dependent NACP savanna plant communities exist within a largely inhospitable matrix of agriculture, plantation forestry, fire-impeding hardwood forest, wetlands, water bodies and urban development, wherein only a small proportion of native savanna species are found (Brudvig & Damschen, 2011;Freeman, Williges, Gardner, & Leone, 2017;Volk et al., 2017). As such, island biogeography concepts are a useful framework for understanding large-scale biodiversity patterns in this region. One such concept is beta diversity, the component of species diversity that accumulates due to site-to-site differences in local species assemblages (Socolar et al., 2016).
Two patterns of beta diversity with profound implications for conservation are turnover (species present at one site are absent at another, but are replaced by other species absent from the first) and nestedness (species present at one site are absent at another, but are not replaced by other species, resulting in smaller assemblages that are nested subsets of larger assemblages) (Baselga, 2010). In landscapes with high species turnover, biodiversity is best conserved by protected area networks that capture the full range of variation within a metacommunity, a conservation strategy that may necessitate protecting smaller and/or more species-poor sites alongside larger and/or richer ones (Socolar et al., 2016). In contrast, networks with high nestedness may benefit from conservation strategies that focus resources on fewer, higher-diversity sites (Socolar et al., 2016).
A related aspect of plant community composition with regional biodiversity implications is the relationship between species richness at small and large scales. Small-scale (≤1 m 2 quadrat) species richness has been a topic of interest among those researching NACP savannas as the first modern plant surveys documented very high plant species richness at these scales (Walker & Peet, 1984). Small-scale richness has been posited as an important conservation target in NACP savannas due to the possibility that species loss at small scales could be a leading indicator of population declines at the whole-preserve level, raising the risk of local extirpation (Glitzenstein et al., 2012;Palmquist, Peet, & Weakley, 2014).
In addition to scaling up to broader spatial perspectives, longer temporal perspectives may also help reveal previously unrecognized patterns and help inform conservation priorities. This is particularly true with respect to the study of fire regime-biota interactions, which are scale-and context-dependent, with outcomes that may only become apparent after many repeated fire cycles (Driscoll et al., 2010;Westgate, MacGregor, Scheele, Driscoll, & Lindenmayer, 2018). Due to the fire-dependent and relatively productive nature of NACP savannas, woody encroachment during the period of intense anthropogenic fire suppression (from the 1930s to the 1990s) was nearly ubiquitous within relict savanna fragments (Heyward, 1939;Van Lear et al., 2005). Managing agencies have since had widespread success in re-establishing the grass-tree vegetation structures thought to be representative of the historical reference condition, through frequent fire and accompanying mechanical and chemical woody reduction treatments. Many savanna preserves throughout the region now experience frequent (<5-year return interval) prescribed fire regimes (Kobziar, Godwin, Taylor, & Watts, 2015), but little is known about the lasting effects of past fire suppression and landscape fragment isolation on the plant species composition of savanna preserves, or how these individual species assemblages fit into the larger regional plant metacommunity.
Like all aspects of fire regime research, understanding the effects of seasonality requires study over long timescales and a variety of contexts. Fire seasonality may have important effects on plant species composition in NACP savannas, as repeated early-and midgrowing season fires result in greater woody mortality and growth reduction (Hmielowski, Robertson, & Platt, 2014), much less vigorous woody resprouting (Drewa, Platt, & Moser, 2002) and higher flowering and seed germination rates among native herbaceous species (Shepherd, Miller, & Thetford, 2011;Streng, Glitzenstein, & Platt, 1993). It has been hypothesized that there may be widespread, undocumented plant adaptations to fires conducted during the months of greatest lightning and wildfire activity (May and June), which is thought by some to constitute the "natural" or pre-anthropogenic K E Y W O R D S beta diversity, fire regimes, functional groups, pine savannas, prescribed fire fire season to which native species would have adapted (Frost, 2006). Fire season variation has also been identified as a potentially important component of conservation strategies, due to phenological differences between species (Hiers, Wyatt, & Mitchell, 2000;Knapp, Estes, & Skinner, 2009), but no studies have examined the plant community effects of NACP savanna fire regimes that are seasonally variable within the same spatial extent.
In this study, we explored fire regime, spatial and vegetation structural drivers of beta diversity and species richness in NACP savannas, focusing on functional groups known to be important for grassland and savanna productivity and biodiversity. NACP savannas are low-N systems, in which some of the most important plant traits with regard to ecosystem functioning are related to N-cycling and availability (Kirkman et al., 2001;Reich et al., 2012). Studies in North American tallgrass prairie communities, which are phylogenetically and functionally similar to NACP savannas, have shown that C 3 graminoids, C 4 graminoids, legumes and non-leguminous forbs exhibit complementarity in root depth distribution and N-acquisition traits, and productivity is maximized when all four of these functional groups are present and species richness within groups is highest (Fargione, Brown, & Tilman, 2003;Fornara & Tilman, 2009;Reich et al., 2012). Woody species are also prominent members of NACP savanna assemblages, which likely add more complexity to the root depth and N-acquisition spectrum, and which have unique direct and indirect relationships with fire regime and understory plant diversity (Veldman, Mattingly, & Brudvig, 2013). We examined functional group richness and beta-diversity patterns within C 3 graminoid, C 4 graminoid, legume, non-leguminous forb (hereafter "forb"), shrub and tree assemblages across a regional sample of 30 different fire-managed savanna preserves in Florida and Georgia, USA. We hypothesized that functional group richness and composition would vary between sites and that some of this variation would be explained by fire regime (return interval, number of fires, time since fire and seasonality), vegetation structure (herbaceous cover, woody cover and tree density) and spatial factors (composition of the surrounding landscape and distance between sites). The scope F I G U R E 1 Locations of the study sites in Florida and Georgia. The subregions (GA south, FL north, FL panhandle and FL central) correspond with differences in mean annual temperature and precipitation due to peninsular climate patterns (PRISM, 2018), which are described further in Appendix 2. The inset map shows the location of the study area within the North American Coastal Plain floristic province of the study also allowed us to test the hypothesis that small-scale richness and whole-site species richness are correlated, both overall and within functional groups.

| Study sites
We sampled the understories of 30 fire-maintained savannas in north Florida and south Georgia, all of which are managed with biological conservation as either a primary or secondary goal ( Figure 1; Appendix 1). Though the longer-term natural and fire history details of any given preserve vary widely, most agencies in the region have been working to implement 2-to 4-year fire return intervals for the past 10-25 years (Kobziar et al., 2015). We sought study sites with little to no history of agriculture other than native range grazing, or forestry practices other than tree removal, as the legacies of more intensive agricultural and forestry activities are known to be associated with novel communities in this region (Brudvig, (Kirkman, Coffey, Mitchell, & Moser, 2004). We also sought sites that had undergone at least three documented fire and/or fire surrogate (mechanical or chemical) treatments, and all but one of the sites met this criterion (the only exception was Osceola National Forest, which was later determined to have had only two recorded treatments). Exotic and weedy species were all but absent from the study sites, which generally typify old-growth savanna communities (sensu Veldman et al., 2015) subjected only to fire suppression.

| Field sampling methods
The study sites were sampled using a nested quadrat sampling design, a modified version of the methods of Peet, Wentworth, and White (1998). At each study site, five randomly placed 10 m × 10 m (100 m 2 ) species-area sampling modules were established. Sampling modules were placed at least 50 m from one another and 20 m from roads or other boundaries. Nested in one corner of each 100 m 2 module were 10, 1 and 0.1 m 2 quadrats. Understory species presence/absence was recorded in all quadrats, and percentage cover to the nearest 5% was recorded in the two smallest quadrats. We classified all understory species by functional group (C 4 graminoid, C 3 graminoid, forb, legume, shrub, tree). C 4 and C 3 graminoid classifications followed Waller and Lewis (1979) and Bruhl & Wilson (2007).
Overstory tree data were collected in the north-west corner of each 100 m 2 module using the point-centred quarter method, from which tree density was then calculated. Field data collection occurred during fall 2014 and fall 2015.
To the greatest extent possible, we located each five-module sampling array within a single burn unit at each preserve, so that all five of the sampling modules shared the same fire and management history. Sampling sites were approximately 9 ha on average, though they varied in shape and size somewhat depending on the shape and size of the surrounding burn unit. At three preserves (Pt. Washington, Bell Ridge and Half Moon), it was necessary to locate sampling modules across more than one burn unit because of size constraints, but in these cases, the units shared the same recorded treatment history.

| Spatial and landscape factors
We stratified the study sites into three soil types, all of which are common NACP savanna substrates, and which are characterized by markedly different plant community compositions (Peet, 2006): Land Use Cover raster dataset (GGISC, 2016). The GLUT land cover categories are much broader than the FLUCCS categories (i.e., evergreen forest and woody wetlands), but were descriptive enough to allow us to roughly calculate the area of each buffer in similar natural vegetation types.

| Analysis of functional group richness
We developed generalized linear mixed models using PROC GLIMMIX in sas v9.4 (SAS Institute, Inc. 2013) to explore the influence of fire regime, spatial and vegetation structure variables on functional group richness. From the fire history questionnaires, we developed a variety of quantitative fuel treatment history parameters, based on histories dating back to 1965 (approximately 50 years from the time of data collection). We began with eleven fire regime variables as potential predictors of functional group richness: total number of fires over 50 years, total number of fire surrogate treatments over 50 years, total number of fire + fire surrogate treatments (weighted equally) over 50 years, mean fire-only intervals over the past 10, 25 and 50 years, mean fire + fire surrogate intervals over the past 10, 25 and 50 years, number of years in active fire rotation and time since fire on the date the data were collected.
We began with five landscape/vegetation variables: number of hectares in same/similar community within 500-and 1,000-m-radius landscape buffers, overstory tree density (trees/ha), percentage shrub cover and percentage herb cover. We included "site" as a random factor in all models and dropped non-significant variables from the models sequentially, in order of least significance, until the lowest AIC was achieved. We calculated the amount of variance explained by each model (R 2 ) using the method of Nakagawa and Schielzeth (2013).
In a separate set of generalized linear mixed models, we analysed interactions and "site" as a random factor. We dropped variables from the models in order of least significance until the lowest AIC was achieved for each functional group and then calculated the total amount of variance explained by each model (Nakagawa & Schielzeth, 2013).

| Analysis of beta diversity
We analysed site-to-site dissimilarity, a measure of beta diversity, between functional group and overall species assemblages in pc-ord version 7 (McCune & Mefford, 2016). We used non-metric multidimensional scaling (NMS) to ordinate the 100 m 2 module data using a Sørenson distance measure (Sørenson, 1948). For the NMS ordination and PerMANOVA, species within each 100 m 2 module were given weighted presence/absence values based on the number of nested quadrat scales at which they occurred (1, 2, 3 or 4). We used PerMANOVA and multivariate pairwise comparisons in PC-ORD, also with a Sørenson distance measure, to test for significant differences in community composition between sites (both overall and within functional groups). We used a Benjamini-Hochberg correction at a false discovery rate of 0.05 to adjust for multiple siteto-site pairwise comparisons.
We performed nestedness analysis using vegan (Oksanen et al., 2013) in R (R Core Team, 2017), with the NODF metric (Almeida-Neto, Guimarães, Guimarães, Loyola, & Ulrich, 2008). We tested the NODF results for significance using the Tuning Peg null model r program (Strona, Ulrich, & Gotelli, 2018), which generates results from 100 null models offering all possible types of row and column constraint, from fully constrained (FF), to partially constrained (PP), to fully unconstrained. FF is the most conservative null model with respect to type I error (falsely detecting nestedness), but has a somewhat higher likelihood of type II error (failing to detect nestedness in a truly nested matrix), while the converse is true for the PP null model (Ulrich & Gotelli, 2012). We report p-values from the fully constrained FF model alongside four incrementally less-constrained null models that lie in the most accurate region of null model space between FF and PP.
We used Betapart (Baselga & Orme, 2012) in R to disentangle the influences of nestedness and turnover on beta diversity, since most real-world metacommunities are characterized by a mix of both patterns (Baselga, 2010). Betapart partitions the Sørenson distance (β SOR , which ranges from 0 to 1) into two additive components, β SIM (Simpson dissimilarity, i.e., turnover) and β SNE (nestedness). These two additive metrics are not direct measures of turnover or nestedness, but rather account for the relative contributions of each phenomenon to the observed overall beta diversity (Baselga, 2010). It is therefore possible for β SNE to be low even when NODF detects significant nestedness, due to a larger turnover component (β SIM ) in the matrix (Baselga, 2010).
β SNE (also denoted β NES in some publications) has an inconsistent and somewhat non-intuitive relationship with absolute nestedness, which is measured most reliably by NODF, and as a result the metric has generated criticism (Almeida-Neto, Frensel, . The biggest weaknesses of the β SNE metric arise when comparing matrices of unequal sizes and in certain limited cases where β SNE can yield false positives (Almeida-Neto et al., 2012;Baselga, 2012). As we were not comparing β SNE between matrices (rather we assessed nestedness and turnover within each functional group × soil matrix separately), and could use NODF results as a benchmark to protect against β SNE false positives, we found the β SOR = β SIM + β SNE calculation provided by Betapart, as well as site-to-site pairwise comparisons (β sor = β sim + β sne ), to be a useful complementary technique. We focused the pairwise analysis primarily on the forb functional group, which is of particular interest because much of the plant diversity of NACP savannas lies in this group and it is among the most vulnerable to species loss via fire suppression (Peet, 2006).
In order to identify drivers of beta diversity, we analysed unfolded Sørenson distance matrices using linear mixed models in the r package lme4 (Bates, Maechler, Bolker, & Walker, 2015). We calculated the Euclidean distance (in km) between sites, as well as the difference between pairs in the following environmental variables determined to be important in the functional group richness models: 25-year fire/fire surrogate total, diversity of seasons, time since fire and 500-m-radius landscape buffer composition. We built linear mixed models to estimate the proportion of variability explained in the Sorenson distances according to the different sources. We included each site in the pair (i.e., site 1 and site 2) as random effects in the models. To estimate the proportion of variation explained we built a null, or intercept-only model (including the random variables), for each soil type. We then built a space, environment and global (i.e., space + environment) model. We calculated pseudo-R 2 based on the proportion of total variance the null model explained compared to the space, environment and global model (Feng, Diehr, Peterson, & McLerran, 2001).
We used Mantel correlograms to quantify the spatial structure in the Sørenson distance matrices, following the geographic lag class method described in Lichstein (2007). We tested for autocorrelation between 100 m 2 modules within different distance classes, the first of which (class 1, <1 km) represents within-site comparisons. The remaining distance lag classes increased by roughly 100-km increments, with some variation in increment size among the different soil types due to the irregular spatial distribution of study sites. Mantel tests were calculated in PC-ORD between each functional group/ soil type Sørenson distance matrix and each binary geographic lag class matrix, yielding a standardized Mantel statistic (r) that was then tested for significance using a permutation test with 2,000 iterations and a Bonferroni correction.

| Analysis of small-scale versus largescale richness
We conducted linear regressions between whole-site species richness and richness at the 0.1, 1, 10 and 100 m 2 levels, in order to determine the degree of correlation between small-scale and wholesite richness.

| Patterns and drivers of functional group richness
We found pronounced differences in functional group richness between sites. Whole-site species richness (the cumulative number of species found across all five sampling modules at each site) varied widely, and the richest sites with respect to any given functional group typically contained anywhere from two to five times more species per functional group than the poorest sites (Table 1).
The final generalized linear mixed model parameters for functional group richness varied by functional group, and there were many significant soil interactions (Table 2) could not be determined from our data. R 2 values were highest for the forb, legume and total models (ranging from 0.66 to 0.75), and lowest in the C 3 graminoid, C 4 graminoid and shrub models (ranging from 0.14 to 0.47).
All three of the seasonality variables (diversity of seasons, growing season fire and lightning season fire) had significant effects, either alone or with soil interactions, in the total richness model (Table 3).
Both diversity of seasons and growing season fire had significant positive influences on forb richness, while the effects of growing season and lightning season fire varied with soil type in the C 3 graminoid and Note. The model for each functional group was selected based on the combination of variables that achieved the lowest AIC. The table shows F statistics and p-values for each included parameter, with significance denoted in bold. The direction of each significant relationship is indicated with a +, − or ± in the case of soil interactions in which the direction differed by soil type. Specifics regarding divergent soil interactions are discussed in the text. The R 2 value is a measure of the percentage of variance explained by each model, derived using the method of Nakagawa and Schielzeth (2013), and the site covariance chi-square is based on a Wald test.
shrub models: positive for Spodosols, negative for Entisols and neutral for Ultisols (Table 3). R 2 values indicate that the fire seasonality models explained roughly half the variation in forb, shrub and overall richness, and very little of the variation in C 3 graminoid richness.
None of the seasonality variables had significant effects on C 4 graminoid richness. Legumes did not have any relationships with seasonal variables, though others have documented sensitivity to fire season in this group (Hiers et al., 2000); this may have been partly due to the fact that our seasonality model contained only Spodosols and Entisols, and the former are relatively inhospitable to legumes.

| Patterns and drivers of beta diversity
Three-dimensional NMS solutions were chosen for all three soil types, with a final stress of 14.0 for the Entisol plot, 13.5 for the Spodosol plot and 12.3 for the Ultisol plot. NMS ordination of the Sørenson distances between overall communities at the 100-m 2 sampling scale revealed that Entisol and Ultisol communities showed distinct separation between sites, even within subregions, while Spodosol sites varied less distinctly (Figures 2-4). Within functional groups, PerMANOVA showed that site was a significant effect for all functional group assemblages on all soil types, and multivariate pairwise t tests showed significant differences between the great majority of site pairs following a Benjamini-Hochberg correction (including within-region pairs), particularly among forb, C 4 grass and legume assemblages (Table 4).
Nestedness analysis using the NODF metric and the Tuning Peg null model algorithm showed that no assemblages were nested at the 0,0 or 1,1 null model positions (the most conservative in identifying nestedness) (Table 5). However, some assemblages tested positive for nestedness from the still-conservative 2,2 position onward (Entisol forbs, Entisol trees and Spodosol C 3 graminoids), and many tested positive at the less-constrained 3,3 and 4,4 positions (shrubs and C 4 graminoids on all soil types), suggesting that some degree of nestedness may be present in those assemblages. There was no evidence of legume nestedness on any soil type. Note. The final model for each functional group was selected based on the combination of variables that achieved the lowest AIC. The table shows F statistics and p-values for each included parameter, with significance denoted in bold*. The direction of each significant relationship is indicated with a +, − or ± in the case of soil interactions in which the direction differed by soil type. Specifics regarding divergent soil interactions are discussed in the text. The R 2 value is a measure of the percentage of variance explained by each model, derived using the method of Nakagawa and Schielzeth (2013), and the site covariance chi-square is based on a Wald test.
Betapart confirmed the Sørenson dissimilarity results derived using PC-ORD, with all functional group × soil type assemblages showing high dissimilarity (i.e., beta diversity) between sites (

| Species richness at different scales
The strongest correlations between quadrat-level and whole-site richness were found among Spodosol sites, particularly in the forb functional group, which had r 2 values of 0.36, 0.48, 0.69 and 0.75 at the 0.1, 1, 10 and 100 m 2 levels, respectively (Table 8). No other soil × functional group combination showed a correlation between whole-site and quadrat richness at the very smallest scale (0.1 m 2 ), and Ultisols also had no correlations between whole-site richness and the next larger scale (1 m 2 ). Entisol forb, legume and total richness at the 1 m 2 level were weakly correlated with whole-site richness.
Correlations between whole-site richness and quadrat richness at the two larger scales (10 and 100 m 2 ) were stronger across the board, particularly at the 100 m 2 level (Table 8). Figure 6 shows a boxplot of 100 m 2 quadrat forb richness superimposed on whole-site forb richness, with sites sorted from highest to lowest mean forb richness at the 100 m 2 level.

| D ISCUSS I ON
We found pronounced variation in functional group richness among NACP savanna sites, with the richest sites typically containing two to five times more species per functional group than the poorest.
Beta-diversity patterns were turnover-dominated for all plant functional groups, a finding that largely held true even among sites in the same climatic subregions of the study area. Though we did find some evidence of nestedness (particularly among Entisol forbs, TA B L E 4 PerMANOVA tests of between-site Sørenson dissimilarity and multivariate pairwise tests of between-site differences in functional group composition (paired t tests with a Benjamini-Hochberg correction) Entisol trees and Spodosol C 3 graminoids), in all cases it was balanced by equivalent or higher turnover, and ultimately, nestedness had only a limited effect on beta diversity for any functional group.
The prevalence of turnover rather than nestedness suggests that even species-poor sites in this portion of the NACP can harbour less-common members of the regional plant metacommunity and may be important targets for conservation despite having smaller assemblages. Our findings indicate that protecting and investing in many sites, regardless of their current species richness, may be the most beneficial conservation strategy for NACP savannas.
The species turnover we observed was due in part to spatial gradients, with many functional group assemblages showing positive autocorrelation in the 1-to 99-km distance band, and most transitioning to negative autocorrelation at 200 km + distances.
Fire regime and landscape composition also explained a portion of between-site variation, though more than half of the observed beta diversity remained unexplained for Entisol and Spodosol assemblages. Some of the unexplained variation likely stems from gradients that were not captured by our set of variables, such as edaphic factors beyond our categorical soil distinctions, or TA B L E 5 Nestedness analysis using the NODF metric, tested against five progressively less-constrained null models generated by the Tuning Peg algorithm Notes. The 0,0 (fixed row totals, fixed column totals) model fully constrains row and column totals, while null models 1,1 through 4,4 allow increasingly more variation in row and column totals, increasing sensitivity to a variety of real-world nestedness patterns but also increasing the risk of type I error. * Significant p-values are denoted in bold.

Nestedness (p-values from five positions in constrained null model space)
TA B L E 6 Betapart partitioning of beta diversity (β SOR ) into nestedness (β SNE ) and turnover (β SIM ) components, broken down by soil type and functional group differences in longer-term fire history beyond the period of record. Still more of the unexplained variation is probably due to long histories of stochastic environmental filtering events not addressed by our study, since droughts, rainfall events and the severity/heterogeneity of individual fires are known to interact with seed arrival to influence NACP savanna community assembly (Myers & Harms, 2011). Finally, a documented lack of dispersal and colonization by native herbaceous species in NACP savannas Kirkman et al., 2004;Ostertag & Robertson, 2007) indicates that these communities are effectively islanded, and the persistence of less-common assemblage members on species-poor fragments suggests a role for neutral immigration and extinction processes in the beta-diversity patterns we observed (Socolar et al., 2016). Our finding that the 500-m-radius landscape buffer composition had a more significant influence on species richness than the 1,000-m-radius buffer composition supports the conclusion that NACP savanna fragments have a limited range of interaction with the surrounding landscape.
As interaction between NACP savanna fragments appears to be relatively limited, the species richness of individual sites may be of critical importance for regional biodiversity. Among Spodosol forb assemblages, we found support for the hypothesis that species richness at small scales is correlated with richness at larger scales, with moderate correlations between whole-site and small-scale (0.1 and 1 m 2 ) richness. On Entisols and Ultisols, however, relationships between whole-site richness and richness at the two smallest scales were weak or absent across all functional groups, a result that we suspect may be due to lower shrub encroachment TA B L E 7 Spatial autocorrelation between 100 m 2 sampling modules at increasing distances, tested using Mantel correlograms at each distance lag class F I G U R E 5 Percentage of beta diversity (as measured by between-site Sørenson distance) explained by null, spatial, environmental and shared spatial-environmental factors, based on linear mixed models. The percentage of variation explained is a pseudo-R 2 based on the proportion of total variance the null model explained compared to the space, environment and global models (Feng et al., 2001) (Glitzenstein et al., 2012). It is important to note that the Ultisol sites in this study were all well drained, whereas other studies in the region have examined plant community change on wetter Ultisol sites, which may be more vulnerable to shrub encroachment (Palmquist et al., 2014). Within-site edaphic heterogeneity is another known correlate of whole-site plant species richness in NACP savannas (Costanza, Moody, & Peet, 2011;Palmquist, Peet, & Mitchell, 2015;Provencher, Litt, & Gordon, 2003) and was likely important on our study sites.
The most consistent fire regime predictors of species richness across functional groups and soil types were the average fire/fire surrogate interval over the past 25 years, and the total number of years a site was in fire rotation. For many variables, only the soil interaction was significant, or the soil interaction was more significant than the main effect. In many cases, the direction of these relationships was the same across soil types, but the slope differed, whereas in several other cases, the relationships varied from positive to negative depending on the soil type. This result was not surprising, as the three soil types differ in many ways that influence vegetation structure, fire regime and the abundance of species within functional groups (i.e., hydrology, edaphic properties, elevations and fuel structures) (Peet, 2006). Of particular note was the positive association between time since fire (up to 4 years) and total Entisol species richness, versus a negative association between time since fire and total species richness on Spodosols and Ultisols, an effect that may point to different woody-herbaceous dynamics on the different soil types. Most of the final models contained current vegetation/landscape parameters as well as fire history parameters, which suggests that the fire regime in previous decades has had effects on species richness that are independent of current vegetation structure, possibly due to within-fragment species extirpation or other compositional changes that occur during periods of fire suppression.
We found empirical support for the theory that fire regimes containing an abundance of growing season fire (March-October) in general, and lightning season fire (May-June) in particular, promote plant species richness in NACP savannas. Diversity of seasons also had a significant positive influence on overall and forb species richness, and these relationships were consistent across soil types.
Though our study did not examine the mechanistic causes of the relationships between burn season diversity and species richness, one logical hypothesis is that burn season variation may better promote co-existence among functional groups and among spring-, summerand fall-flowering forb species. It may be that prescribed fire regimes with an emphasis on lightning season fire, but also containing varia- The continued application of fire in the coming decades will likely foster conditions hospitable to higher species richness, but evidence from this study and others suggests limited dispersal of native species between fire-maintained fragments Kirkman et al., 2004;Ostertag & Robertson, 2007

APPENDIX 4
Means and 95% confidence intervals for 100 m 2 species richness parameters used in generalized linear mixed models