Novel insights into how the mean and heterogeneity of abiotic conditions together shape forb species richness patterns in the Allegheny plateau ecoregion

Abstract While plant community theory tends to emphasize the importance of abiotic heterogeneity along niche axes, much empirical work seeks to characterize the influence of the absolute magnitude of key abiotic variables on diversity. Both magnitude (as reflected, e.g., by a mean) and heterogeneity (variance) in abiotic conditions likely contribute to biodiversity patterns in plant communities, but given the large number of putative abiotic drivers and the fact that each may vary at different spatiotemporal scales, the challenge of linking observed biotic patterns with the underlying environment remains acute. Using monitoring data from a natural resource agency, we compared how well statistical models of the mean, heterogeneity, and both the mean and heterogeneity combined of 17 abiotic factor variables explained patterns of forb species richness in Northeast Ohio, USA. We performed our analyses at two spatial scales, repeated in spring and summer across four forest types. Although all models explained a great deal of the variance in species richness, models including both the mean and heterogeneity of different abiotic factors together outperformed models including either the mean or the heterogeneity of abiotic factors alone. Variability in forb species richness was mostly due to changes in mean calcium levels regardless of forest type. After accounting for forest type, we were able to attribute variation in forb species richness to changes in the heterogeneity of different abiotic factors as well. Our results suggest that multiple mechanisms act simultaneously according to different aspects of the abiotic environment to structure forb communities, and this underscores the importance of considering both the magnitude of and heterogeneity in multiple abiotic factors when looking for links between the abiotic environment and plant community patterns. Finally, we identify novel patterns across spatial scales, forest types, and seasons that can guide future research in this vein. OPEN RESEARCH BADGES This article has earned an Open Data Badge for making publicly available the digitally‐shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.5061/dryad.kp3cb17.

However, it is often difficult to incorporate realistic abiotic complexity into research questions because organisms can respond to multiple abiotic factors in different ways. For example, species in a community might show classic niche differentiation with respect to one abiotic variable (McKane et al., 2002), while simultaneously sharing a single niche optimum for some other variable (Vergnon, Ooi, & Freckleton, 2017). Individual abiotic factors vary through time and have unique spatial characteristics, differing in their range of variation and scale of spatial autocorrelation. Given all this complexity, it is easy to see why the role of abiotic drivers is often simplified, especially in research emphasizing other aspects of plant community structure. Establishing how plant communities are impacted by multiple abiotic factors that vary at different spatial and temporal scales is nevertheless a crucial step toward a more complete picture of spatial plant ecology.
Recent simulation work has shown promise in using spatial statistics to distinguish between plant community patterns that arise from different mechanisms (Brown, Illian, & Burslem, 2016). However, as the authors point out, ignoring the complexity of the abiotic environment will likely limit the application of this approach to real systems, because the signature that results from a given biotic mechanism may depend on the spatial structure of the abiotic environment. In other words, regardless of whether our central interest is in the abiotic drivers themselves or in other processes that structure plant communities, a comprehensive understanding of the influence of abiotic conditions is critical.
In these situations, abiotic conditions at a site could be measured in terms of the mean level of each important factor, and the species richness at that site results in part from how many species have these mean conditions within their niche (Scheffer & van Nes, 2006).
If the total number of species able to tolerate different abiotic conditions varies in a systematic way (Myers & Harms, 2009), such that there is a tendency toward greater richness as conditions become more favorable, then we would expect to see a correlation between the mean abiotic conditions and species richness. Relationships between mean conditions and richness may be further influenced by relative fitness differences among competitors for different levels of a limiting resource (Clark & Tilman, 2008), or by interactions with dispersal ability leading to mass effects (Mouquet & Loreau, 2003;Shmida & Wilson, 1985).
Abiotic heterogeneity-as measured by the spatial variance or coefficient of variation of abiotic variables-can also be related to patterns of species richness (Stein, Gerstner, & Kreft, 2014).
Relationships between species diversity and abiotic heterogeneity are commonly referred to as heterogeneity-diversity relationships (HDRs; Lundholm, 2009), terminology that we adopt throughout this paper. We might expect plant richness to increase with abiotic heterogeneity if communities are structured by niche differences (Chu & Adler, 2015;Kneitel & Chase, 2004;McKane et al., 2002).
HDRs can also arise in other ways, such as when abiotic heterogeneity promotes microbial diversity which in turn increases plant species diversity (van der Heijden et al., 1998), or if sampling units encompass a wide range of abiotic conditions, for example, across habitat types with different regional species pools (Davies et al., 2005).
There is no a priori reason that both the AE and HDR hypotheses cannot act simultaneously (HilleRisLambers et al., 2012;Kerr & Packer, 1997;Kreft & Jetz, 2007). When they do, both the mean and the variance in abiotic conditions will influence plant diversity.
Observational studies can be extremely helpful for confronting the complexity inherent in the sheer number of distinct ways that abiotic conditions can scale up to community patterns. Using observational data from a natural resource agency, we show how characterizing both the mean and heterogeneity in abiotic conditions simultaneously is a powerful way to identify potential mechanisms structuring herbaceous plants across spatial scales, seasons, and forest community types. This type of analysis can be a productive precursor to future experiments that pin down mechanism, particularly in understudied systems.
The current study addresses the following three gaps common in plant community ecology: (1) plants respond to an array of abiotic factors (Pausas & Austin, 2001), but many studies only include one or a few factors (Stevens & Carson, 2002).
(2) When studies do include more than one abiotic factor, they are usually concerned with the effect of either abiotic mean conditions or abiotic heterogeneity on community patterns (e.g., Burton et al., 2011;Tamme, Hiiesalu, Laanisto, Szava-Kovats, & Pärtel, 2010; but see Richard, Bernhardt, & Bell, 2000), or they treat mean conditions and heterogeneity as alternative hypotheses (Lundholm, 2009). And finally, (3) Abiotic drivers change across spatial and temporal scales (e.g., soil calcium with bedrock type [Bellemare, Motzkin, & Foster, 2005], or nitrogen with season, as in the vernal dam hypothesis [Tessier & Raynal, 2003]), but many studies consider only a single scale (Baer, Blair, Collins, & Knapp, 2004;Kreft & Jetz, 2007). We addressed points (1)-(2) through a statistical analysis of 29 herbaceous layer forb communities located across four common forest types found in the Allegheny plateau ecoregion, for which we have estimates of both the mean and spatial variance of 17 different abiotic factors (34 abiotic factor measurements). To explore point (3), we ran analyses across 2 spatial scales, repeated in spring and summer.
Using these data, we built a suite of linear and linear mixed-effects models to determine which abiotic factor measurements explained the most variation in species richness within and across different forest community types at two spatial scales in spring compared to summer.

| Study system
Forbs-flowering, nongrass herbaceous plants-are an important yet understudied component of forest ecosystems (Gilliam, 2014; Figure 1). Often both common and diverse, forbs contribute to ecosystem services like nutrient cycling (Gilliam, 2007), and serve as nectaries and hosts for pollinators (Bess, 2005;Davis & Cipollini, 2014;Hanula, Ulyshen, & Horn, 2016). They also account for a large proportion of plant biodiversity in temperate forests. For example, most spring ephemerals (species that begin to bloom in early April, then die back and give way to a summer vegetation community at the end of May) are forbs, and therefore drive much of the compositional change in plant communities across seasons (Gilliam, 2014).

| Study area
Sampling was conducted in collaboration with a monitoring program implemented by Cleveland Metroparks (CMP). The sites used in our analysis are a subset of 100 sights chosen randomly using the Generalized Random Tessellation Stratified (GRTS) spatial sampling method at the start of the monitoring program in 2010 (Mack & Robison, 2010). At each site, a 400 m 2 plot composed of four contiguous 100 m 2 subplots (of any configuration) was placed to encompass an area with similar canopy species, and without any conspicuous inclusions (e.g., a vernal pool, or a drastic change in slope aspect).
One plot with just two subplots was also included in our analysis.
Spring sampling was conducted by a crew led by the first author (17 April-14 May), and summer sampling was conducted by two crews of CMP employees led by the second author (23 June-8 September). CMP's crews only sample in summer, so our spring sampling was done without the benefit of CMP's field staff. For this reason, visiting all 100 sites was not feasible in spring. As such, we selected 4-5 plots in each of 6 land holdings (29 plots total) that had at least one tree individual >60 diameter at breast height (DBH; to exclude early successional communities). In each land holding, we attempted to select plots representing the four terrestrial hardwood community types commonly found in Northeast Ohio: Beech-Maple, Floodplain, Mixed, and Oak forests (Eysenbach & Hausman, 2013), although this was not always possible ( Figure 2).

| Field sampling protocol
Sampling procedures were based on the North Carolina Vegetation Survey Method (Peet, Wentworth, & White, 1998). Data were collected on all plant species (including non-natives) found growing to the height of the tallest herbaceous species off the forest floor, including tree seedlings, ferns, grasses, sedges, mosses, and forbs. As we did not directly control for functional and/or phylogenetic relatedness in our study, we limited our analyses to include only forb species because we were especially interested in spring ephemeral communities. Species cover was assessed in each subplot (α) and each plot (γ) using the Braun-Blanquet method (Peet et al., 1998;Figure 2). Cover data were converted to presence-absence data to get a species count in each plot (plot-level species richness, γ S ) and subplot (subplot-level species richness, α S ). Sometimes, individuals of the same genus but different species (identifiable via unique characteristics, like smooth or hairy stems in the genus Solidago) were nevertheless unidentifiable to species at the time of sampling. These were identified as Genus sp. 1, 2, 3, etc., allowing us to accurately measure species richness without a confirmed species identification.
Approximately 340 g of soil taken from within 10 cm of the surface were collected with a soil probe at randomly chosen points throughout each subplot (8 plugs/subplot), providing subplot-level composites for each abiotic factor j, α j . Undried soil samples were stored in Whirl-Pak bags at room temperature for the sampling season, and then sent to A&L Great Lakes Laboratories (Fort Wayne, Indiana, USA) where samples were dried, ground, sieved and then analyzed for a suite of variables known to affect plant performance: organic matter (%), phosphorus (Bray 1 parts-per-million [ppm]), pH, cation exchange capacity (meq/100 g), nitrogen (% total), carbon (%), and the carbon:nitrogen ratio. Potassium, magnesium, and calcium were also measured, as both absolute quantities (e.g., in ppm), and relative to one another (e.g., % base saturation; Table 1).
In the center of each subplot, leaf litter and organic layer depth were measured to the nearest 10th of a centimeter, and light data were collected using a densiometer. Restrictive layer depth was also measured in the center of each subplot using a tile probe (Table 1).
F I G U R E 1 Example forb community in a forest understory located in Northeast Ohio. Species pictured include Cardamine concatenata, Dicentra sp., Geranium maculatum, Floerkea proserpinacoides, Podophyllum peltatum, and Viola sp.
Organic layer depth in summer was extremely zero-inflated and was excluded from the analysis.

| Seasonal designations
Each plant was designated a spring or summer species based on the date of up to three flowering specimens collected since 1990 found digitized in the Ohio State Herbarium (Museum of Biological Diversity Herbarium, The Ohio State University, 2016; Table 2).
Observed flowering times were also recorded in the field by the first author and used to corroborate and supplement dates collected online. Spring and summer species had an average flowering date before and after June 1, respectively (except Allium tricoccum, which was designated as a spring species since it loses its spring vegetative parts before flowering and senescing in June [Vasseur & Gagnon, 1994]). Only designated spring species were used in the spring sample, whereas all identifiable plants were used in the summer sample, and thus often included spring ephemeral species that had not yet senesced, especially early in the summer season. We chose this approach because spring species that persist into summer have fully emerged, while summer species that start to appear in spring have not. Species that were only identifiable to genus were included in the analysis and designated a season based on the natural history of their taxonomic group (e.g., violets only identifiable to Viola sp. were designated as spring species).

| Community patterns and abiotic measurements
Species presence data were recorded for each plot (to compute plot-level species richness, γ S ) and subplot (to compute subplotlevel species richness, α S ; Figure 2). As all abiotic factors were measured at the subplot level, we were able to assess the magnitude of abiotic conditions at both the subplot level and, by averaging across subplot-level values, at the plot level (α j and ᾱ j , respectively). Abiotic heterogeneity could only be assessed at the plot level and was measured as the coefficient of variation (CV) across subplot-level values within a plot (cv(α j )) ( Figure 2).
Abiotic factors were measured in units of very different magnitude (Table 1). To facilitate comparison between linear models, all abiotic predictor variables (i.e., the raw data taken at the subplotlevel, the plot-level means as measured via the raw data, and the plot-level coefficient of variation as measured via the raw data) were first normalized, if strongly skewed, via a log or square root transformation. All measurements were then standardized (ztransformed) to have mean of 0 and a variance of 1.

F I G U R E 2
The distribution of sampled plots and the forest community types they represent within Cleveland Metropark's park boundaries (located around Cleveland, Ohio). Pictured alongside is an example of a single plot, with levels γ (whole plot) and α (subplots) used to denote at what scale measurements were taken 2.6 | Analysis

| Statistical modeling and model selection
Our first goal was to identify the best statistical models of plot-level and subplot-level species richness (γ S and α S , respectively) via model selection. Our general approach was to begin with a full model for species richness, γ S , that included the plot-level means of all 17 abiotic factors (ᾱ j ) and then conduct model selection to determine the best predictors of species richness. Except where noted, all model selection was done using the dredge function in the R package MuMIn (Bartoń, 2018), wherein subsets of a full model were ranked according to their AIC C score (Hurvich & Tsai, 1989). Since there was collinearity between multiple abiotic measurements ( Figure 3), we restricted the dredge so that any pairs of abiotic terms with >40% correlation would not be included together in any model subsets.
After each dredge, models with an AIC C score within 2.0 of the best model's AIC C were retained. Diagnostics were then run on TA B L E 1 Minimum and maximum measured values of all predictor variables α j for abiotic factor j across all subplots α

Spring Summer
Organic matter (%) α OM 1.7-14.6 1.7-31.5 Phosphorus (  We next repeated this procedure with a full model that included the plot-level heterogeneities (cv(α j )) instead. Once the best meanbased and CV-based models were attained, we combined the surviving terms from each into a third full model and conducted model selection again, to see whether the combined model would better predict species richness compared to the best models based on abiotic mean or CV alone. To facilitate comparison of AIC C scores, we removed the same outliers from all full models within a season (e.g., if a plot was an outlier in the mean model, it was also removed from the CV model). We identified the best statistical models for subplotlevel species richness (α S ) in the same way, using sampled subplotlevel conditions (α j ) as abiotic predictors.

| Assessing different model structures
The procedure described in the preceding paragraph was repeated for several classes of "full model," including linear models that treated all plots as independent sampling units, and linear mixed-effects models that allowed plot or subplots in different forest types to differ in their slope and/or intercept (Table 3). In other words, linear mixed-effects models were used to account for the nonindependence of plots found in similar forest communities (step B in Table 3).
Each class of model allowed us to answer a different question about the richness patterns we observed. Linear models (step A in Table 3) describe the relationship between species richness and abiotic factors across all plots, regardless of forest community type.
Linear mixed-effects (LME) models (steps C-D in Table 3) describe the relationship between species richness and abiotic factors within forest types. For example, a variable intercepts model (step C1) suggests that mean richness across forests varies but the overall relationship with the abiotic factor does not. A variable slopes model (step C2) suggests that the overall relationship with an abiotic factor depends on forest community type. In steps D1 and D2, best models from C1 and C2 were modified to allow slopes and intercepts to vary simultaneously.
We did not build linear models for subplot-level richness because of the spatial proximity of subplots within plots. Instead, we compared two types of LME models: one with subplots grouped by plots, and another with subplots grouped by plots nested in forest type.
Models were screened for variable slopes in an iterative fashion (i.e., separate models were built to cycle through each abiotic measurement separately, instead of dredging a full model). In models with variable slopes and intercepts (e.g., step D1 in Table 3), we did not force slopes and intercepts to vary according to a particular abiotic factor. This allowed us to see whether a different abiotic measurement was important within different forest community types, compared to those driving patterns across plots and/or forest types.
To make sure this practice was not obscuring a best model, we also built expanded models where any abiotic measurement associated with variable slopes was also added to the suite of best predictors identified in previous steps. In all LME models, a diagonal positive definite matrix was used to specify the variable slopes part of the model to control for any covariance between intercepts and slopes.
All linear models were fit using ordinary least squares. All LME models were fit using maximum likelihood in the lme function in the R package nlme (Pinheiro, Bates, DebRoy, Sarkar, & R Core Team, 2018).
The variation explained by LME models is given in terms of marginaland conditional-r 2 values (Johnson, 2014). The marginal-r 2 value is the variance explained in species richness across plots-for example, the single intercept and slope for each predictor that describes its effect on species richness without accounting for forest type. The conditional-r 2 value in the LME model is the marginal-r 2 plus the variance explained by the intercepts and/or slopes that vary by forest type.
Each step in Table 3 yielded best candidate models describing how richness is related to mean abiotic conditions, abiotic heterogeneity, or both (except B, which considered the null model including no abiotic predictor variables). From among these, AIC C scores and the amount of variance explained across study levels were used to choose a suite of overall best models for each season. All analyses were conducted in R (R Foundation for Statistical Computing).

| Variance partitioning
After assembling a suite of best models, variance partitioning was used to determine whether more variance in plot-level species richness could be attributed to mean abiotic conditions or to abiotic heterogeneity. This allows a higher-level synthesis of the results from our series of fitted models. For each season, the varpart function in the R package vegan (Oksanen et al., 2018) was used to partition the variance in species richness explained into the portion attributed to mean abiotic conditions, the portion attributed to the heterogeneity in abiotic conditions, the portion attributed to both (i.e., indistinguishable due to correlation), and the portion that remained unexplained by our abiotic measurements.

TA B L E 2 (Continued)
Using the abiotic predictors from the subset of best model(s) from the statistical analysis, we did variance partitioning at two different levels: once using all plots and excluding abiotic predictors whose slopes varied with forest community type, and once for plots within each forest type and including abiotic factor predictors whose slopes varied depending on forest. In the latter analysis, we only partitioned variance between two abiotic predictors at a time since sample size was low within forests (n = 6, 9, 5, and 9 for Beech-Maple, Floodplain, Mixed, and Oak forests, respectively). We then compared results across all pairwise combinations of abiotic predictors, but retained only those pairs that had the smallest amount of unexplained variance to use in our analysis.

| Statistical modeling
Mean abiotic conditions and heterogeneity in abiotic conditions across multiple factors together explained substantial variation in plot-level forb richness in both seasons (Table 4). For example, across plots (i.e., from step A in Table 3), plot-level richness increased with mean calcium availability in both seasons (Figure 4a

| Variance partitioning
Across all plots, most of the variability in plot-level species richness was attributed to mean abiotic conditions in spring, and in summer this trend became even more pronounced in that no variability could be attributed to abiotic heterogeneity alone (Figure 6a). Within most forest community types, however, more variability in species richness was attributed to abiotic heterogeneity than to mean abiotic conditions, and this was true for both seasons (Table 5, Figure 6b,c).

| D ISCUSS I ON
In spring, 50% of the variation in forb species richness was explained by forest community type alone, with no measured abiotic predictors included in the linear mixed-effects models. In summer, only 30% of the variation in forb species richness was explained by forest community type, likely because Floodplain and Mixed forests, as well as Oak and Beech-Maple forests, had similarly sized species pools in summer. In other words, mean richness varied in a predictable way across forest community types in spring, and somewhat less so in summer. Differences in the mean number of species could be due to differences in the size of a shared species pool (e.g., if certain forest types include additional species in an otherwise shared pool), or to the presence of a unique species' pools for each forest community type. Notably, the AICc score of more complicated models including abiotic predictors were always much lower than the AICc scores of the null model including only forest community type (∆AIC C = 15.9 in spring and ∆AIC C = 16.2 in summer). A lower AICc score compared to the null model indicates that the additional variation in forb species richness patterns explained after adding abiotic predictors, for example, the intercept(s) and slope(s) of the models shown in Table 4 are not simply due to a more complex model, but actually reflect additional variation explained by the measured abiotic factors after accounting for differences in mean richness across forest types. We review which predictors best explained forb species richness and offer some brief interpretation of our results based on when we found support for the AE and HDR hypotheses.

| Summary
Our results confirm that patterns of forb richness depend on both mean abiotic conditions and abiotic heterogeneity (Table 4, Figure 4).
In spring, including both the mean and heterogeneity of different abiotic factors resulted in the best overall model and explained the most variation in plot-level richness (Table 4), giving simultaneous support for the presence of heterogeneity-diversity relationships (HDRs) and the available energy (AE) hypothesis.
In summer, the mean model was the most supported according to its AIC C score, but the heterogeneity model and the combined model TA B L E 3 Stepwise protocol used to generate candidate best models for species richness Step

Model type Question Modeling plot-level richness (γ S )
A: across plots Linear Can abiotic factors explain species richness across plots? explained more variance (Table 4). We suspect that the heterogeneity and combined models have inflated AIC C scores because in summer, the penalty for dividing forest types into 4 categories was not justified given that Mixed and Floodplain forests had the same mean species richness, as did Oak and Beech-Maple forests. Since there was little correlation between means and CVs of the abiotic measurements included in the summer models (Figure 3), we think the explanatory power of the combined model is biologically meaningful and not simply due to collinearity. We therefore interpret our results as evidence that summer plot-level richness is also determined by both mean abiotic conditions and their heterogeneity.
Because our analysis makes use of an observational dataset, we cannot isolate the mechanisms underlying our results. Nevertheless, our findings are suggestive of which mechanisms are likely to be structuring community patterns at different study levels. Below, we highlight a few intriguing results and discuss possible mechanisms that might explain them.

| Positive heterogeneity-diversity relationships: pH and nitrogen heterogeneity in summer
In summer, richness increased with increasing pH and nitrogen heterogeneity: both positive HDRs (Figure 4g,h). This suggests that assemblages of summer species may differentiate along pH and nitrogen niche axes, while the lack of a similar HDR in spring suggests that assemblages of spring species do not. It is also possible that the spatial structure in pH and nitrogen changed seasonally such that HDRs occurred at a different scale in spring than summer, allowing only the summer relationship to be detected by our sampling scheme.
Many herbaceous forbs found in these forests form and benefit from associations with arbuscular mycorrhizal fungi (AMF; Burke et al., 2018;Whigham, 2004). Soil chemistry properties including pH and nitrogen heterogeneity are likely to influence mycorrhizal communities themselves (Coughlan, Dalpé, Lapointe, & Piché, 2000;Dumbrell, Nelson, Helgason, Dytham, & Fitter, 2010;Kluber et al., 2012), as well as the relationship between plants and mycorrhizae (Johnson, Graham, & Smith, 1997). More research is needed to assess whether plant-microbe interactions could scale up to contribute to the positive HDRs we observed in summer.

| Negative heterogeneity-diversity relationships: potassium and phosphorus heterogeneity in spring
In spring, we found that plot-level richness decreased with increasing potassium and phosphorus heterogeneity: both negative HDRs (Figure 4e,f). Negative HDRs can arise due to dispersal limitation (Lundholm, 2009). If spring species that differentiate across potassium and phosphorus conditions are also dispersal limited, an inability to arrive and fill available niches could lead to low richness in plots with high abiotic heterogeneity. Alternatively, if more heterogeneous plots contain subplots that can support fewer species, this too would explain the TA B L E 4 Selected models explaining plot-and subplot-level species richness |forest " indicates what grouping factor was used in linear mixed-effects models-that is, either forest, plot, or plots in forests (i.e., " |forest∕plot "). The value listed in the r 2 column is either the standard multiple-or adjusted-r 2 in linear models with one predictor or multiple predictors, respectively, or the marginal-r 2 for linear mixed-effects models.
F I G U R E 4 Selected trends in plot-level richness (γ S ). Plotted relationships correspond to models listed in Table 4 (and described in Table 3)
appearance of negative HDRs (Yang et al., 2015). Indeed, mean phosphorus in spring is positively correlated with phosphorus heterogeneity ( Figure 3). This means heterogeneous plots tend to include higher phosphorus concentrations, which could conceivably speed competitive exclusion (Tilman, 1982) and lead to lower plot-level richness (e.g., in Figure 4c,d,j).

| The available energy hypothesis: calcium, nitrogen, and other measures of habitat quality
When species respond to the magnitude rather than the heterogeneity of an abiotic factor, increasing the magnitude or "available energy" (AE) of that factor may either increase or decrease species richness. Theory supports both of these possibilities: increasing the availability of a limiting nutrient can lead to faster competitive exclusion (Clark & Tilman, 2008;Tilman, 1982), as mentioned above.
Alternately, increasing the level of a factor that is favorable to many different species, can increase the diversity of species that can thrive at that location (Lundholm, 2009). Our results cannot speak to mechanism, but we observed both patterns. Richness decreased with increased nitrogen availability in beech and Oak forests in spring (Figure 4j), and across plots in summer with increasing potassium and phosphorus availability (Figure 4c,d). Conversely, in spring, plot-level richness increased across mixed forests with increasing nitrogen (Figure 4j). Forb richness also increased with increasing calcium availability in both seasons and at both spatial scales ( Figures   4a,b and 5a,b). The strength and persistence of this particular trend throughout our study, as well as its prevalence in the literature (Bellemare et al., 2005;Burton et al., 2011;McEwan & Muller, 2011;Peet et al., 2014), suggests it is worth further consideration. We note that mean calcium is also highly correlated with pH and potassium variability ( Figure 3).
F I G U R E 5 Selected trends in subplot-level species richness (α S ). Plotted relationships correspond to models listed in Table 4. Numbers indicate to which plot each subplot belonged. The x-axes for all abiotic predictors have been normalized for comparison, and community pattern response variables that were log-or square root-transformed for analysis have been plotted with the original values on the y-axis F I G U R E 6 Proportion of variance explained in plot-level species richness, γ S , partitioned between mean abiotic conditions (light gray), abiotic heterogeneity (CV, in dark blue), both (mean + CV in white), and neither (unexplained in black). To partition variance explained across all plots before accounting for forest type (a), we excluded abiotic factor predictors that changed trends across forest community types (e.g., abiotic factors in column 3 of Figure 4). To partition variance explained within forest community types (b and c), we partitioned the variance explained in plot-level species richness in each forest and included abiotic factor predictors that changed trends across forest community type

| Effects of scale
To determine whether the response to resource availability depended on spatial scale, we also identified the best statistical models for subplot-level species richness using sampled subplot conditions (α j ) as abiotic predictors and asked how results changed compared to plot-level relationships. In both seasons, subplot-level richness and plot identity were strongly correlated (step B conditional-r 2 = 90% in spring and summer). Plot-and subplot-level richness were also strongly positively correlated (r 2 = 84% and 88%). In other words, species-rich plots had species rich subplots. Since species richness was so strongly correlated at the plot-and subplot-levels, and since calcium availability was an important predictor of plot-level richness, we thought that calcium availability would also be related to forb richness at the subplot level.
Interestingly, we found that changes in mean calcium within plots did explain additional variation in species richness (compared to models without abiotic predictors in step B; ∆AIC C = 12.0 in spring and 24.7 in summer); however, this relationship was much stronger in summer (slope = 0.73) compared to spring (slope = 0.37; Figure 5a,b). In spring, subplot-level richness did increase with increasing calcium (Figure 5a), but the relationship was weak, and most of the subplot-level variance was only explained by accounting for plot identity (marginal-r 2 = 11%, conditional-r 2 = 88%; Table 4). Thus, in spring, spatial scale mattered: calcium availability explained patterns of plot-level richness better than patterns of subplot-level richness. In summer, both plot-and subplot-level richness increased with increasing mean calcium (Figures 4b   and 5b). What could cause seasonal differences in scale-dependent effects of calcium is an interesting question for future studies.

| Variance partitioning and the importance of forest community type
In both seasons, the variance in species richness across plots was attributed largely to mean abiotic conditions, regardless of forest type ( Figure 6a). The variance in species richness explained within each forest, however, was attributed to the heterogeneity in abiotic conditions of different factors (Figure 6b,c and Table 5).
In spring, most of the variability in species richness was attributed to potassium heterogeneity in both Beech-Maple and Floodplain forests (Figure 6b and Table 5). In summer, most of the variability in species richness was attributed to pH heterogeneity in Beech-Maple and Oak forests, and to heterogeneity in cation exchange capacity in Mixed forests (Figure 6c and Table 5). Thus, we find evidence that plot-level abiotic heterogeneity can drive patterns of species richness (Figure 6b,c), but that the identity of the important abiotic factor may change with season and forest community type (Table 5).

| Limitations and future directions
Patterns of richness and composition within forest type could be confounded by differences in size, location (e.g., more urban or rural), and/or invasive species and deer population density (Eysenbach & Hausman, 2013) across land holdings. To control these confounding variables and test proposed mechanisms, it is necessary to confirm the drivers of increased species richness experimentally.
We also cannot discern from our data whether species are stably coexisting (due, e.g., to niche differentiation) or whether they are simply co-occurring (due to mass effects, high available energy, slow competitive exclusion, etc.). Pending further empirical investigation, we think it is useful to consider coexistence and species richness, although clearly related, as potentially arising via separate processes.
There was a great deal of correlation between abiotic factors in both spring and summer (Figure 3), which we addressed by allowing only uncorrelated (<40%) combinations of abiotic factors into a model, thus excluding information about potentially confounding correlations. As such, there were many nearly equivalent models during the model selection process. Qualitatively, our results provide a framework for understanding the joint contributions of abiotic conditions across plots and abiotic heterogeneity within forests to patterns of forb species richness ( Figure 6), but larger sample sizes within forest type are needed, as are experiments that disentangle correlated variables (e.g., between calcium availability and pH, calcium availability and potassium heterogeneity; Figure 3). To clarify how species richness is structured by different abiotic factors, we also hope that the relationships found at different study levels will be able to guide future experiments designed to test the hypotheses presented above.

| CON CLUS IONS
The abiotic component of what shapes herbaceous layer forb communities is complicated, but the large amount of variation explained across study levels justifies further study and discussion of this complexity. We underscore the importance of considering mean abiotic conditions and heterogeneity simultaneously across multiple abiotic factors. In particular, we show that patterns of TA B L E 5 Pairs of abiotic measurements within each forest type that resulted in the lowest amount of unexplained variation in plotlevel species richness after partitioning variance due to the abiotic mean, the abiotic CV, both, or neither  species richness vary across forest types according to abiotic quality, as well as within forest type according to abiotic heterogeneity.
Connecting patterns associated with the abiotic mean and heterogeneity back to theoretical insights allowed us to suggest which mechanisms could be structuring plant communities at different spatial scales and seasons. Future empirical work in herbaceous forb communities should focus on using experiments to pin down the mechanism underlying plot-level increases in richness with increasing nutrient availability, and investigating potential interactions between biotic processes and abiotic drivers within plots; for example, whether abiotically mediated microbial communities could be driving fine-scale patterns of forb richness in different ways in different forests.

ACK N OWLED G M ENTS
The authors would like to acknowledge the seasonal field crews at Cleveland Metroparks for collecting the summer data, and to Andrew Lance for his help collecting the spring data. We would also like to thank members of the Abbott Lab, plant ecologists at Cleveland Metroparks, the Case Western Reserve University Ecology & Evolution Journal Club, as well as Robert K. Peet and three anonymous reviewers for invaluable input on previous drafts.

CO N FLI C T O F I NTE R E S T S
The authors report no conflict of interests.

AUTH O R CO NTR I B UTI O N S
S. Catella conducted the spring survey and S. Eysenbach conducted the summer survey. All statistical analyses were drafted and conducted by S. Catella. All authors contributed to the drafting and revising of the final manuscript.