Simulating Groundcover Community Assembly in a Frequently Burned Ecosystem Using a Simple Neutral Model.

Fire is a keystone process that drives patterns of biodiversity globally. In frequently burned fire-dependent ecosystems, surface fire regimes allow for the coexistence of high plant diversity at fine scales even where soils are uniform. The mechanisms on how fire impacts groundcover community dynamics are, however, poorly understood. Because fire can act as a stochastic agent of mortality, we hypothesized that a neutral mechanism might be responsible for maintaining plant diversity. We used the demographic parameters of the unified neutral theory of biodiversity (UNTB) as a foundation to model groundcover species richness, using a southeastern US pine woodland as an example. We followed the fate of over 7,000 individuals of 123 plant species for 4 years and two prescribed burns in frequently burned Pinus palustris sites in northwest FL, USA. Using these empirical data and UNTB-based assumptions, we developed two parsimonious autonomous agent models, which were distinct by spatially explicit and implicit local recruitment processes. Using a parameter sensitivity test, we examined how empirical estimates, input species frequency distributions, and community size affected output species richness. We found that dispersal limitation was the most influential parameter, followed by mortality and birth, and that these parameters varied based on scale of the frequency distributions. Overall, these nominal parameters were useful for simulating fine-scale groundcover communities, although further empirical analysis of richness patterns, particularly related to fine-scale burn severity, is needed. This modeling framework can be utilized to examine our premise that localized groundcover assemblages are neutral communities at high fire frequencies, as well as to examine the extent to which niche-based dynamics determine community dynamics when fire frequency is altered.

Introduction

Wildland fire is a globally critical process for maintaining biodiversity in many terrestrial ecosystems (Pausas and Ribeiro, 2017). Surface fire regimes are particularly important, where frequent low-intensity fires often maintain the highest levels of floral and faunal diversity (Bond and Keeley 2005; Mitchell et al., 2009; Pausas and Ribeiro, 2017). These frequently burned ecosystems are globally widespread and of significant conservation concern (Staver et al., 2011); thus, understanding how fire structures patterns of diversity is a critical need. In particular, longleaf pine (Pinus palustris Mill.) woodlands of the southeastern US coastal plain are dependent on frequent fire with burns occurring as frequently as every 18 months. This ecosystem is a reflection of the fire regime, where frequent fires support both survival and regeneration of the dominant overstory tree, the longleaf pine, and keep other vegetation predominantly in understory stature (Boyer, 1990; Jose et al., 2006). Across the range of longleaf, prescribed burning is applied at 1- to 5-year return intervals to maintain the structure and function of these endangered ecosystems (Jose et al., 2006; Mitchell et al., 2006).

A unique characteristic of these systems is the high plant species richness and endemism within the groundcover plant community that occurs at fine scales and is linked to fire frequency (Hardin and White, 1989; Barnett, 1999; Dell et al., 2019). The relationship between high fire frequency and high groundcover diversity is well known (Peet and Allard, 1993; Taylor, 1998) and historically documented (Small, 1923, Wahlenberg, 1946), but the mechanisms driving this connection still remain elusive (Kirkman and Myers, 2018). Previous work (Platt et al., 1988; Kirkman et al., 2001; Pausas and Ribeiro, 2017) has focused on documenting correlations between species diversity, fire return intervals, soil edaphic gradients (i.e., productivity), and seasons of burn. Other studies have tested the competitive effects of dominant guilds (bunchgrasses or shrubs) on rare species in these frequently burned systems and found no competitive effects (Myers and Harms, 2009), while others found facilitation effects (Iacona et al., 2012). These studies do not, however, account for microscale spatial demographics that might structure the groundcover community and do not attempt to project richness patterns based on demographic parameters. Furthermore, many small plant species are found in high densities on very uniform sandy soils (up to 40–50 within 1 m2, Peet and Allard, 1993; Kirkman et al., 2001). If niche-based processes were structuring these communities, resource-partitioning for dozens of species would occur within just a few meters across these resource-poor soils.

Research on ecological community assembly and biodiversity has been invigorated by the unified neutral theory of biodiversity (UNTB, Hubbell, 2001; Chave, 2004; Matthews and Whittaker, 2014; Missa et al., 2016). The theory has stimulated modeling and empirical studies examining the relative roles of stochastic processes versus competitive exclusion and niche differentiation in driving patterns of species richness and community assembly, particularly among plants (McGill, 2003; Volkov et al., 2003; Chave, 2004; Wootton, 2005; John et al., 2007; Rosindell et al., 2011). UNTB assumes that individuals of different species within a trophic level are functionally equivalent and posits that communities of these species are structured by random demographic processes, such as birth, death, and dispersal limitation, and at longer timescales, speciation, and extinction.

For many ecosystems, both neutral processes and niche differentiation affect community structure (e.g., Gravel et al., 2006; Adler et al., 2007; Matthews and Whittaker, 2014), with the UNTB treated primarily as a null model that captures the impact of unknown mechanisms (Alonso et al., 2006; Rosindell et al., 2012). We hypothesize that the high species richness of frequently burned longleaf pine ecosystems, and potentially other similar systems, may be an actual consequence of these neutral processes at smaller temporal and spatial scales because fire plays two crucial roles. (1) Repetitive fire (1- to 5-year intervals) suppresses competition of fast-growing woody plant species (Barnett, 1999; McGuire et al., 2001) and acts as a competitive filter (Pausas and Verdú, 2008; Myers and Harms, 2011), and (2) fire causes random mortality at fine scales (Wiggers et al., 2013; O’Brien et al., 2016a). These random mortality events are a defining process of the UNTB. Fire-driven mortality in these systems occurs at the scale of an individual plant (Wiggers et al., 2013; O’Brien et al., 2016a), because the relevant scale of variability in fire intensity of these surface fires is also found at this same fine scale (Hiers et al., 2009; Loudermilk et al., 2012). These neutral dynamics, however, are usurped when fire is removed from these systems as predictable changes in the plant community occur, with faster growing fire-sensitive species released and competition drives the shift to woody plant dominance (Brockway and Lewis, 1997; Glitzenstein et al., 2003). Ultimately, the relative importance of mortality and other neutral parameters (birth and dispersal limitation) in these or similar groundcover communities associated with fire adapted forests remains unexplored.

Despite criticisms and attempts to reject the UNTB, the utility of estimating and examining variation in demographic parameters in the context of neutral contributions to biodiversity is clear (Rosindell et al., 2012). Whether or not it is treated as a null model, UNTB allows for deviations from neutrality, and different demographic parameters are more sensitive to change. Thus, it is important to examine how changes in neutral parameters affect diversity predictions across different communities and across different scales (Hubbell, 2001; Matthews and Whittaker, 2014). Properly quantifying variation in neutral parameters within and among modeled communities and at multiple scales requires careful development to ensure accurate predictions of species richness. There are a number of approaches for examining such variation, including sensitivity analyses, which may be used to quantify the relative differences in sensitivity between parameters (Hamby, 1994; Saltelli et al., 1999). These analyses are used with newly developed models to provide insight into how parameters vary and to generate guidance on appropriate model interpretation and future directions.

In this study, we explored how neutral processes (Hubbell, 2001) might be the mechanism maintaining high groundcover diversity when exposed to frequent fire. We explored UNTB-based parameters within a modeling environment using empirical observations. We utilized 4 years of plant census data with individual plants mapped at very fine scales (100 cm2), following their fates before and after two low-intensity prescribed burns in a northwest FL, USA, longleaf pine ecosystem. We estimated UNTB-based parameters for two parsimonious autonomous agent models, which were distinct by spatially explicit and implicit local recruitment processes. We employed parameter sensitivity tests within our two models to examine how parameter estimates affected output species richness. We also examined how varying input species frequency distributions and community size assumptions of the model influenced parameter sensitivity across models.

Materials and Methods

Study Site

This study was conducted at Eglin Air Force Base (EAFB), FL, from 2012 to 2016. EAFB, the former Choctawhatchee National Forest, is located in northwest FL, USA (30°34′23.50″N, 86°34′48.64″W). The site includes nearly 180,000 ha of longleaf pine forests, which contains over half of the remaining old growth stands of this forest type (Varner et al., 2000; Holliday, 2001). The study sites were within the Southern Pine Hills District of the Coastal Plain Physiographic Province with deep, well-drained sandy soils (Overing et al., 1995). The climate of the area is subtropical with mean annual temperatures of 19.7°C and mean annual precipitation of 1,580 mm (Guo et al., 2016). Elevations were 52 to 85 m above sea level. Vegetation was dominated by a longleaf pine overstory, with a sparse midstory of various deciduous oaks, for example, Quercus laevis Walter, Quercus margaretta Ashe, Quercus incana Bartram, and Quercus germinata Small.

Data Collection and Prescribed Burning

We established our plots within existing EAFB monitoring plots in two habitat types—sandhills and flatwoods. The sandhills habitat is xeric with deeper sands and at a higher elevation than the mesic flatwoods habitat. The flatwoods have more productive soils than the sandhills. In each habitat type, we established 15 random plots measuring 1 m × 3 m for a total of 30 plots distributed across nine management units. Every fall from 2012 to 2016, yearly records of all living plants within each plot (sandhills n = 15 and flatwoods n = 15) were mapped using a grid coordinate system measuring 10 cm × 10 cm. All plant individuals were recorded by species. The resulting plant dataset includes four intervals of data comprising 123 plant species and over 7,000 individuals. A list of plant species is found in . Experimental prescribed burns were conducted in Spring 2013 and 2015 across all plots. The fires were conducted as part of normal forest management operations and burn prescription conditions at EAFB (O’Brien et al., 2016a; O’Brien et al., 2016b).

Autonomous Agent Modeling

To groundcover community dynamics, we developed one spatially explicit model and one spatially implicit neutral autonomous agent simulation model, respectively, the “spatial model” and the “non-spatial model,” using the python programming language (v. 2.7.11 Python Software Foundation). These two models were distinct strictly by their recruitment processes; they were designed to assess the degree to which spatial and non-spatial dispersal functions influenced simulated species richness patterns. As the input data were collected in frequently burned areas of EAFB with a 2-year return interval, we assumed that fires were implicit in both models. This is realistic because after 2 years without fire, species richness begins to decline in these ecosystems (Brockway and Lewis, 1997; Glitzenstein et al., 2003). A conceptual model of model inputs and processes is found in .

Figure 1

Conceptual diagram of this study’s neutral model used to simulate groundcover species richness in a frequently burned landscape. The metacommunity is the groundcover species frequency distribution within each longleaf pine habitat (flatwoods and sandhills) as recorded in this study at Eglin Air Force Base (EAFB). These data were used as input at five scales (described in text) to initialize each community simulated in the model. Each individual experiences a probability of mortality. Then empty cells experience a probability of birth, where species are recruited from either the local community pool or metacommunity pool, with the latter through immigration. The box inset describes how local recruitment is simulated in the spatial and non-spatial versions of the model. The resulting communities are used as the starting communities for the subsequent time steps. *The size of each community and number of communities are determined by the community size parameters: area and number of areas. ^If immigration occurs during a “birth” event, the metacommunity pool is used in place of the local community (community) pool for recruitment.

UNTB requires a metacommunity species pool, a local community species pool, and community level birth, death, and immigration (dispersal limitation) rates. To estimate the distribution of plant species found across each habitat, we combined all plant data from the burned areas to create one input species distribution for each habitat (). These metacommunity species distributions (“the metacommunity”) were used to populate the initial distribution of plant species within each simulated area by randomly choosing species from this distribution. The metacommunity also provided a source for recruitment that occurred from immigration only (not the local community) throughout each simulation.

We simulated individual plants through a presence-absence approach, where each individual was subjected to a yearly mortality risk and empty cells to recruitment potential. The probability of birth (recruitment) and death were produced using a random number generator and a set seed with no density dependence effects on death. We used a lattice-based modeling approach with a 1-year time step, where only one individual could occupy a cell in the lattice. Each cell was 100 cm2, the same scale at which individuals were sampled in the field. The spatial model was distinguished from the non-spatial model by the implementation of periodic boundaries, which eliminated edge effects on spatially dependent recruitment.

For the non-spatial model, local recruitment was a function of the probability of birth, where species were sampled from the frequency distributions for each simulated area (i.e., the local community pool). For the spatial model, local recruitment was a function of the probability of birth and the local neighborhood of individuals. Here, we used a negative exponential function within the NumPy (v1.11.3) module for python, to determine which individual from the neighboring community will be recruited into the focal cell. We set a maximum distance of five cells in all directions from the focal cell for local dispersal to standardize the local recruitment pool. This sets a maximum dispersal distance to 0.71 m, the distance diagonal from the focal cell. This includes an area up to 1.1025 m2 (1.05 m × 1.05 m), which represents the local community pool at any given “birth” event. Although seed dispersal is quite variable among these groundcover species (Stamp and Lucas, 1990), we found this seeding distance to be comparable with that of a similar longleaf pine reference site, where most species had restricted seed dispersal mechanisms (Kirkman et al., 2004).

Model Parameterization and Sensitivity Analysis

Mortality and birth rates were estimated from the plant monitoring data. We used plant data from the 2012 to 2014 and 2014 to 2016 surveys, because its period was a consistent 1.5-year time period between each burn and the following survey (e.g., Spring burn 2013 and Fall 2014 survey), including the 2012 survey, which initially burned in Spring 2011. An individual plant of an identified species was considered dead if it was not found in its original cell in the next census. Recruitment occurred when a plant was newly found in a cell in the next census. Mortality rates (1|0) and birth rates (1|0) were averaged across plots, years, and habitats. Sensitivity analyses were then utilized because these resulting mortality and birth rates were highly variable, and when these values were applied to the models, they had confounding effects on model outputs. We employed the Fourier amplitude sensitivity test (FAST) as defined by Saltelli et al., (1999) to provide an unbiased approach to assess the effects of each range of parameter values using a concise method. The FAST was applied in the python library SALib (v.1.0.2) to examine the first-order effects of model parameters on simulated species richness. FAST quantifies the sensitivity of multiple parameters simultaneously by assigning each input parameter a frequency at which to oscillate across a range of input values. Decomposing the model outputs using these same frequencies allows for the extraction of the influence of each input parameter on the model results. This allows model sensitivity to be assessed using far fewer simulations than a complete exploration of the parameter space. By using this method, a unique set of input parameter values was created for each model run using the empirical data.

Five model parameters were used in the FAST; the three UNTB-based or neutral parameters (mortality, birth, and immigration) and two more exogenous-type parameters that influence the spatial boundaries in which the model runs, namely, area (size of each community) and number of areas (or communities). The two additional parameters (i.e., community size parameters) were included because they were found to drive community size during model initialization, recruitment, and immigration, all influencing output species richness. Input parameter sets were created for each of these five model parameters within each habitat.

The values used as input parameter sets in the FAST include birth rate (0.2–0.6), mortality rate (0.2–0.6), immigration rate (0.01–0.25), area simulated (1–25 m2), and number of areas simulated (1–30 areas). We chose a range of plot (area) sizes and numbers of plots (areas) used to estimate species richness in these small plant communities (Hardin and White, 1989; Peet and Allard, 1993; Kirkman et al., 2001; Dell et al., 2017) to explore their effects on simulated outputs. Our immigration rates are a range of values that represent a high migration system (>0.1) and low migration system (<0.1), which affects rates of species turnover (Volkov et al., 2003). In SALib, 100 samples for each of the five parameter sets were created. This resulted in 500 unique parameter combinations used to simulate each habitat and scale of input for the spatial and non-spatial models. To include variability, three replicates using three different random number seeds were simulated. Our target output variable was “normalized richness,” the proportion richness change from time step 1 to the ending time step of 50, when model surpassed any instability.

Furthermore, we explored issues associated with initial conditions to assess their impact on simulated species richness. This was performed because it has been found that initial conditions, particularly related to scale of the species frequency distributions, can influence simulations of neutrality (Hubbell, 2001; Matthews and Whittaker, 2014). As such, we examined effects of the scale of inputs by coarsening the initial species frequency distributions using a presence–absence approach across various resolutions, resulting in inputs at the scales of 10 cm × 10 cm (original resolution), 20 cm × 20 cm, 50 cm × 50 cm, 100 cm × 100 cm, and 200 cm × 200 cm. Within each resulting coarser size cell, species presence was recorded. These species distributions () were used as initial inputs to assess effects of scaled initial conditions on the FAST simulations.

In SALib, the simulation outputs were analyzed across the 500 parameter combinations and 20 scenarios (2 models × 2 habitats × 5 scales of inputs). We focused on the first-order effects of the parameters, illustrating the percent of the variance in normalized richness, which is explained by each parameter alone. This allowed us to focus on individual parameter sensitivity to a multi-faceted combination of initial conditions and model processes.

Results

Our autonomous agent models simulated groundcover species richness distributions across two longleaf pine habitats ( and ). The results revealed there was little difference in outputs between habitats ( and , and ), which is a result of the similarity in size and shape of their input species frequency distributions or their metacommunities (). For brevity, we focus on the sandhills habitat results here but include all flatwoods outputs in and . The FAST illustrated that the direct first-order effects of varying mortality and birth parameters (from 0.2 to 0.6) contributed to less than 5% of the model variance (). Immigration was the most sensitive of the neutral parameters, varying from about 5% up to 20%. The community size parameters (number of areas and area) were similarly sensitive, varying also up to 20%, and illustrated the most sensitivity between replicates (). There were little differences between species richness distributions and parameter sensitivities between the spatial and non-spatial models ( and ).

Figure 2

Normalized species richness of the autonomous agent model at various scales. Normalized species richness is the species richness change from time 0 compared with the last time step (year 50) of each model run using the sandhills habitat data. Normalized species richness is illustrated between scales of input data and spatial versus non-spatial dispersal models.

Figure 3

First-order effects (%) of the three neutral parameters (mortality, birth, and immigration) and spatial boundary parameters (area and number of areas) on simulated normalized richness using different scales of input data from the sandhills habitat using the spatial model (A) and the non-spatial model (B). The legend illustrates the resolution (in dm) at which the data were created from the empirical data, including the original resolution (1.0: 10 cm × 10 cm). Boxplots represent variance across three replicates of the model run using 500 FAST parameter combinations.

The model was sensitive to the scale of the input frequency distributions () and had mixed effects on individual parameter sensitivity (). At the finest scales of input distributions (10 cm × 10 cm), immigration rate (compared with birth and mortality) had the largest individual parameter sensitivity on normalized richness (18–20% vs. <3%). Generally, as scale of input distributions coarsened, the neutral parameter sensitivities decreased, particularly for immigration. Effects of scale on mortality rates were minor, with little variability across scales of input data and low overall individual parameter sensitivity. Scale had significant consequences for the sensitivity of both community size parameters. In general, as the input scale coarsened, the individual parameter sensitivity increased, rather than decreased, as seen with the neutral parameters. The FAST also illustrated some nonlinear differences in sensitivity across scales for the community size parameters (e.g., “area” parameter in ).

Discussion

This modeling experiment provided support for neutral processes as a mechanism driving groundcover community dynamics in this high fire frequency ecosystem. This is the first study to examine explicitly fine-scale plant demographics; prior research focused on correlations at the habitat-level (Kirkman et al., 2001; Pausas and Ribeiro, 2017) and plot-level manipulation experiments on species competitive traits (Myers and Harms, 2009; Iacona et al., 2012). These findings point to the importance of neutral processes in supporting the high species diversity in longleaf pine ecosystems, and they help explain the tight coupling among fire frequency, species richness, and the abatement of competition. This study also supports the central assumption of the UNTB: that niche-based processes do not drive community assembly in these groundcover communities when fire frequency is high.

We simulated a range of groundcover species richness distributions () using simple probability functions, representative of the UNTB. We found that a broad range of estimates for mortality and birth (values from 0.2 to 0.6) can be used to simulate species richness. This suggests that robust estimates are produced without including individual species and their respective traits, and efforts to obtain better estimates may be unnecessary. However, the high sensitivity of dispersal limitation (immigration) illustrates the importance of quantifying accurate immigration rates and of understanding community level dispersal when modeling fine-scale community richness. The importance of immigration, effects on species turnover, the difficulty in its parameterization, and the underlying assumptions were similar to those of other UNTB simulated systems (Chesson, 2000; Girdler and Barrie, 2008; Lowe and McPeek, 2014; Matthews and Whittaker, 2014) and should be explored further.

We found that simulated community stability and rarity were influenced by the community size parameters that affected community size and feedbacks with the neutral parameters (). Simulating greater number of areas or larger spatial extents yielded more stable communities and a greater likelihood for rarer species to persist through time when communities are large. The number of areas simulated directly affected the number of local communities, while area simulated affected size of the local community pool. Furthermore, the low effects from different functions of dispersal (spatial vs. non-spatial dispersal, and ) suggest that defining dispersal processes may be trivial when modeling species richness at these fine scales. In reality, dispersal distance is important for recruitment processes that generate and maintain understory plant diversity. As such, restricting recruitment in the local community using the spatial model can create a modeling environment conducive to realistic dispersal patterns of these small plants (Stamp and Lucas, 1990; Kirkman et al., 2004). Here, various dispersal distances and functions can be manipulated within the model and explored across sites. These spatial dispersal processes may be particularly important at coarser scales, where dispersal distance is more of a relevant process, which applies to mosaics of forest management and fire regimes across the landscape.

Community structures, represented here as different input species frequency distributions derived from different scales, were also important ( and ). At coarser scales, the input frequency data were less skewed, where rare species were no longer comparatively rare and simulated richness was more stable than using finer scale data. The effects of scale on each individual parameter were complex (). At the finest scales of input distributions, immigration rate had the largest individual sensitivity among the neutral parameters affecting normalized richness. Here, rarer species depended more on immigration from the metacommunity than using coarser scaled inputs. The lower individual sensitivities of birth rate illustrated that rarer species were more dependent on effective recruitment at finer scales than at coarser scales. The overall low sensitivity of mortality illustrated that at even high mortality rates (here, up to 0.6), the community would thrive if the local and metacommunity pool was abundant and available. As such, even if birth rates were low, high mortality values would provide more space—open cells—available for new recruits.

This study illustrated the relative importance of our neutral model parameters and how community structure and spatial boundaries can influence simulated species richness. Our next step is to critically examine the multi-scalar outputs to the empirical data and examine which parameter combinations best predict high species richness. This study focused on first-order effects, yet second-order effects, or inter-parameter sensitivity can also be explored further to examine the more nuanced and complex effects of parameter interactions on model outputs. Furthermore, these results do not illustrate the spatial or temporal scales at which neutral processes might break down, specifically where there are changes in fire frequency, land use, overstory structure, fuel characteristics, and soil properties, as well as variability in fire intensity within and among burns. These effects of scale on community structure and rarity are, however, consistent with both neutral theory and island biogeography predictions (Volkov et al., 2003; MacArthur and Wilson, 2016) and provide an example of patterns of diversity maintained by both processes. Furthermore, the fine-scale quantification of frequency distributions for these groundcover species represents the spatial scale (100 cm2) at which community dynamics occur, which was important for examining model simulations and parameter sensitivity.

An alternative approach of characterizing over 100 plant species by their life history traits and incorporating them into a niche-based model (Peterson, 2006) would be significantly more complex than this simple modeling approach. This would likely provide little to no improvement in predicting species richness patterns under high fire frequencies, particularly if stochastic processes are driving community dynamics. Utilizing species traits would be useful, however, when quantifying the extent to which niche-based processes determine community dynamics when fire frequency is altered and in examining the impact of variation of fire intensity and resulting burn severity within and among burns (O’Brien et al., 2016a; O’Brien et al., 2016b). A deeper understanding of the impact of variable fire regimes would require combining a neutral model with a competitive exclusion model operating at the functional group level. This strategy has been explored when examining the impact of stochastic versus deterministic processes across ecosystems (defined as the Continuum Hypothesis, Gravel et al., 2006). For this more inclusive framework, deterministic and neutral assumptions form the ends of a gradient, defined in this case by fire frequency or spatial scale. However, we argue that the continuum hypothesis plays out within a single landscape driven by the dynamics of fire. When fire is removed, the transition from neutral to deterministic community dynamics can also vary across landscapes, with competitive exclusion occurring more quickly in higher productivity sites than on less productive xeric sites (Kirkman et al., 2001).

Frequently burned ecosystems represent a reservoir of biodiversity worldwide (Noss et al., 2015), and it is critical to understand the processes by which fire maintains diversity. This will aid in supporting critical elements of managed fire regimes, setting realistic restoration targets, and examining monitored changes in plant communities (Hiers et al., 2012). These microscale groundcover community dynamics could be important through their interaction with larger scale disturbances, such as vehicle-driven soil compaction, machine-constructed fire breaks, and timber management practices. The extent and spatial pattern of each disturbance will likely influence dispersal patterns that ultimately affect predictions of establishment and assumptions of neutrality (Gardner and Engelhardt, 2008; Lowe and McPeek, 2014). For example, timber harvesting causes changes in understory fuels, which could have feedbacks on how the understory responds (changes in fire mortality) to the subsequent altered patterns of fine-scale fire intensity (Mitchell et al., 2009). We found that simple neutral-based simulations were useful for projecting species richness at multiple scales, although further work is needed to determine optimal parameterization, including dispersal variation, and to test against empirical data. Examining model transferability across fire-prone ecosystems is an important next step, particularly because this neutral model is scalable; the model is dependent on ecosystem species frequency distributions, empirical neutral parameters, and implicit fire effects on richness. Furthermore, the model was intentionally constructed using cellular automata properties that can be integrated with similarly constructed models of forest fuels (Hiers et al., 2009) and fire behavior (Loudermilk et al., 2012; Achtemeier, 2013). In fact, our overarching goal is to create a spatially explicitly unified model of fire behavior and fire effects, which includes simulating fire impacts on groundcover community assembly. Ultimately, these results are globally relevant because similar woodland and grassland ecosystems that are structured by the keystone process of frequent fire occur worldwide, and all are threatened by alterations to their fire regime. Examples include tallgrass prairies of North America (Leach and Givnish, 1996), African savannas (Bond et al., 2005), Australian eucalypt woodlands (Bowman, 2000), and cerrado savannas of Brazil (Mistry, 1998). These fire-maintained ecosystems may be governed by neutral processes generated by fine-scale fire mortality and fire-driven competitive abatement to maintain diversity in their groundcover communities.

Data Availability

All datasets generated for this study are included in the manuscript and the , and the raw plant species data at https://github.com/spokswinski/groundcoversim.

Author Contributions

EL, LD, JO’B, SP, AH, BH, LR, JD, SG, and JH contributed to manuscript writing and advised in the development and interpretation of the autonomous agent model, sampling design, and analyses. EL was the lead in model development, FAST analysis, and manuscript writing, with significant contributions by LD and JO’B. All authors contributed to idea development and implementation. Other contributions include, but are not limited to, field plant monitoring by SP; sensitivity analysis by SG; and prescribed burning by BH and JH.

Funding

This research was funded by the US Department of Defense Strategic Environmental Research and Development Program (#RC-2243) and the National Science Foundation (DEB-1442103). This work was also supported by the US Department of Agriculture Forest Service National Fire Plan.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.