Temporal stability of grassland metacommunities is regulated more by community functional traits than species diversity

Temporal stability in aboveground net primary productivity (ANPP) is influenced by several attributes of plant communities. Identifying the primary regulators of stability and their roles across spatial scales is of both practical and theoretical importance. We assessed effects of species diversity (local or alpha diversity and species dissimilarity between local communities) together with spatial differences in two community functional attributes (mean aboveground biomass and community leaf dry matter content, LDMC) on temporal stability in spring ANPP of restored grassland. Biomass, community LDMC, and species dissimilarity were derived from remote measurements of canopy reflectance of grassland on two soil types. Results demonstrated that productivity at the larger spatial scale of the metacommunity (communities connected by dispersal) was stabilized more by spatial differences in community functional traits than by diversity or community differences in diversity. Spatial differences in community biomass and LDMC stabilized metacommunity productivity by increasing differences in the productivity responses of spatially distinct communities to interannual variation in precipitation, but de‐stabilized ANPP on one soil type by reducing the temporal stability of local communities. Our results demonstrate the utility of remote sensing for quantifying community attributes useful to assess or predict temporal stability of grassland ANPP. We conclude that temporal stability in productivity depended largely on community differences in functional attributes that couple plant growth to changes in resource availability.


INTRODUCTION
Research into the influence of biodiversity on ecosystem temporal stability has intensified in the wake of concerns over habitat degradation and environmental changes. Understanding how various components of biodiversity influence stability and how their influences may vary across spatial scales is of obvious relevance to fostering sustainable ecosystems. Much of biodiversity-stability research has focused on relationships between plant species diversity and temporal stability in biomass production, where stability typically is defined as the ratio of mean biomass production to its temporal standard deviation (µ/σ). Heightened species diversity increases biomass temporal stability via a variety of nonexclusive mechanisms that increase mean biomass production, reduce the response of biomass to perturbations, or some combination (Loreau et al. 2001, Isbell et al. 2009, Hector et al. 2010, Gross et al. 2014, Tilman et al. 2014). It has become clear, however, that the biodiversity effect is complex, often mediated through ❖ www.esajournals.org 1 July 2020 ❖ Volume 11(7) ❖ Article e03178 multiple mechanisms , Craven et al. 2018. Species diversity alone may account for relatively little of variation in temporal stability in some ecosystems (S. Wang, M. Loreau, C. de Mazancourt et al., unpublished manuscript). Several potentially correlated aspects of biodiversity contribute to the full suite of statistical and biological mechanisms that promote temporal stability. Craven et al. (2018) identified four classes of biological drivers of the biodiversity-stability relationship: species richness, functional diversity, phylogenetic diversity, and functional composition. We focus on two of these drivers, plant species diversity and plant community functional composition as reflected in community-weighted functional attributes. Both species diversity and community functional attributes affect temporal stability through pathways that are linked to species differences in trait values. Species diversity includes variation in values of all traits and trait combinations in a community that are coupled to taxonomic diversity. Community attributes, by contrast, depend primarily on trait values of dominant species, consistent with the mass ratio hypothesis (Grime 1998).
Trait differences underpin effects of both diversity and community attributes on temporal stability, but the two biological drivers of the biodiversity-stability relationship need not be strongly correlated. Garnier et al. (2004), for example, reported a negative correlation between specific aboveground net primary productivity (SANPP) and community-weighted values of leaf dry matter content (LDMC) among succession communities, but found no correlation between SANPP and the number of species that accounted for 80% of maximum standing biomass. The implication is that community LDMC and species richness were not correlated among communities. In general, we expect a stronger correlation between species diversity and community attributes when diversity differences arise along disturbance gradients or gradients in abiotic variables that simultaneously influence species composition and the values of favored traits. Severe disturbances, for example, favor both low species diversity and convergence of species into similar functional groups (Biswas and Mallik 2011).
Most manipulative studies of temporal stability have focused on local processes. Efforts increasingly are being directed at understanding regulation of temporal stability at larger spatial scales, including the scale of the metacommunity. As defined here, metacommunities are aggregate communities consisting of two or more communities interconnected by dispersal. The development of a partitioning model of temporal stability has significantly aided the quest to better understand stability regulation across spatial scales (Fig. 1). Loreau (2014, 2016) partitioned temporal stability of ecosystems at large scales into two multiplicative components. Metacommunity temporal stability (γ stability; γ S ) is the product of temporal stability of component communities (alpha stability; α S ) and asynchrony among communities (asynchrony; ω). Asynchrony in productivity occurs when plant growth responds differently to drivers of production among spatially distinct communities. The stability partitioning approach follows from the partitioning model of regional diversity into alpha diversity (α D ) and spatial beta diversity (β D ; Whittaker 1972). Each component of γ S potentially is influenced by both community attributes and diversity.
Asynchrony among communities, like γ S , usually is assumed to vary with species diversity. Asynchrony typically is analyzed as a function spatial variation in plant species composition (β diversity) on the premise that differences in species identities or abundances promote differences in the response of plant production to environmental fluctuations. Increasing spatial variation in species composition can increase temporal stability in productivity (McGranahan et al. 2016, Wilcox et al. 2017S. Wang, M. Loreau, C. de Mazancourt et al., unpublished manuscript) as well as enhance other aspects of ecosystem functioning (Pasari et al. 2013, Van Der Plas et al. 2016, Hautier et al. 2018. However, there is evidence that the asynchrony effect on γ S depends, sometimes primarily, on factors in addition to species diversity. Spatial asynchrony in productivity should be maximally responsive to spatial differences in those factors that most directly link growth dynamics to the primary drivers of ANPP variability (Loreau and de Mazancourt 2013). At the spatial scale of tens of hectares, these growth-regulating factors may better be ❖ www.esajournals.org 2 July 2020 ❖ Volume 11(7) ❖ Article e03178 reflected in differences in the functional composition of communities than in diversity per se. McGranahan et al. (2016), for example, found that temporal variability at larger spatial scales in tallgrass prairie declined (γ S increased) as synchrony in dynamics among functional groups of species decreased (asynchrony increased). Similarly, metacommunity stability depended more on asynchronous fluctuations in populations of the same species among communities than on β diversity (Wilcox et al. 2017). Diversity and trait-based community attributes are providing insight into the regulation of temporal stability across spatial scales but remain a challenge to measure frequently and over large areas. Airborne remote sensing techniques for plant biochemical and biophysical features are increasing capacity to quantify stability and its vegetation correlates at expanded spatial scales. Remote sensing techniques now are routinely applied to estimate biomass production as required to calculate temporal stability, for example. Techniques also are being developed to estimate community attributes that may explain local and spatial variation in ecosystem stability. Remote measurements of vegetation have, for example, been used to estimate communityweighted values of plant traits linked to ecosystem function (Chadwick and Asner 2016, Ali et al. 2017, Polley et al. 2020) and components of species diversity (Dalmayne et al. 2013, Möckel et al. 2016, Polley et al. 2019b).
Here, we assess relationships between spatial components of the temporal stability of ANPP in restored grassland (α S and spatial asynchrony) and two classes of vegetation properties. Properties considered include species diversity (alpha diversity, α D , and species dissimilarity between communities as an index of β D ) and two community attributes (mean biomass production and community LDMC). Community LDMC is a functional attribute (trait) that mechanistically links variation in plant species abundances and ecosystem functioning (Duru et al. 2013). Species dissimilarity between communities, α D , and the abundance-weighted LDMC of grassland were derived from hyperspectral measurements of canopy reflectance (Polley et al. 2019b(Polley et al. , 2020. Precipitation variability is the primary driver of interannual variability in ANPP in the grassland studied (Polley et al. 2020).
Our specific objective was to assess contributions of two spatial components of diversity (α D and dissimilarity) and two community attributes (biomass production and community LDMC) or attribute differences between communities to temporal stability of metacommunity ANPP during spring (γ S ) using grassland communities Fig. 1. We assessed the contributions of community attributes and species diversity to spatial components of the temporal stability of metacommunity aboveground net primary productivity (γ S ). Spatial components of γ S include temporal stability at the local scale (alpha stability of communities in the metacommunity; α S ) and spatial asynchrony between communities (asynchrony; ω). Variation in α S and ω was analyzed as a function of spatial differences in two community attributes, mean aboveground biomass and community leaf dry matter content (LDMC), and two diversity components, either alpha diversity (α D ) or species dissimilarity between communities.
❖ www.esajournals.org 3 July 2020 ❖ Volume 11(7) ❖ Article e03178 growing on the same or different soil types. Metacommunity stability may be jointly controlled by diversity and community attributes. For example, ANPP stability of communities of perennial species depended on both diversity and a community attribute (weighted SLA; Polley and Wilsey 2018). Diversity and community attributes, in turn, can depend on soil properties (Polley et al. 2019b(Polley et al. , 2020. We used structural equation models to assess interactive effects of diversity and community attributes on temporal stability of metacommunities on the same or different soil types. Little is known of the contribution of community attributes to regulation of temporal stability across spatial scales. This study addresses this knowledge gap. Consistent with previous findings, we predicted that (1) variation in γ S and its spatial components would depend as greatly on community attributes as on diversity (Polley et al. 2013, Polley andWilsey 2018), and (2) linkages between stability and both diversity and community attributes would differ between soil types.

Site
We measured temporal stability in aboveground productivity during spring over 4 yr in restored grassland located in Temple, Texas, USA (31°10 0 N, 97°34 0 W; 223 m above sea level). Annual precipitation averages 880 mm (92 yr record), with peaks in spring and autumn. Maximum values of mean monthly temperature vary between 15.3°C and 35.4°C in January and August, respectively.
Eight randomly selected stands, each 17-m wide and 137-218 m long (0.26-0.37 hectares), were restored to grassland in 2010 by planting a mixture of 38 native perennial forb and grass species (equal weight of seed per species; Appendix S1: Table S1) in a former agricultural field ( Fig. 2; Long-Term Biomass Experiment; LTBE). Stands traverse a catena from a silty clay soil on the upland (Austin series; 43% clay) to a clay soil on the lowland (Houston Black series; 52% clay). Soil types are distinguishable in space. Remaining stands (16) in LTBE were planted to a monoculture of the C 4 grass Panicum virgatum L. (cultivar "Alamo"). Thirteen 7 m diameter patches were permanently located along the upland to lowland catena in each stand of restored grassland, resulting in a total of 52 patches per soil type. Restored grassland is not grazed or fertilized but is hayed annually following the growing season. Ungrazed grassland in the area typically is hayed annually.
The composition of restored grassland is strongly influenced by a cohort of volunteer annual grass and forb species that invaded following restoration. Most of these species establish following autumn rainfall and complete growth during spring. Dominant species of annuals include the grass Bromus japonicus Thunb. ex Murray and forbs Gaillardia pulchella Foug and Monarda citriodora Cerv. ex Lag.

Field measurements
We estimated cover by species and of bare soil in eight randomly chosen grassland patches on each soil type during both April-May and August 2017 (32 patches total). Cover per patch was calculated from visual estimates of species and bare soil cover in eight randomly located plots (76-cm diameter) in each patch (Fig. 2). Species diversity per patch was calculated as the exponential of Shannon entropy (Shannon diversity) for i = 1 to S species, where p i = proportion of total plant cover for species i across all plots per patch. The exponential transformation converts Shannon entropy to its number equivalent, reflective of true diversity (Jost 2007). Biomass per unit of surface area for each patch was estimated by averaging aboveground harvests from two plots per patch. Relationships between temporal stability and both species diversity and community functional attributes were investigated using metacommunities created by aggregating grassland patches (Fig. 2). We created 48 communities from the 104, 7-m patches measured. Communities were created by pairing adjacent patches in each stand (n = 13 patches per stand), beginning with the first two patches on the upland of each stand and proceeding from upland to lowland with patches 3 and 4 through 12 and 13. One patch per stand was omitted prior to creating communities in order to equalize the total number of ❖ www.esajournals.org 4 July 2020 ❖ Volume 11(7) ❖ Article e03178 Fig. 2. We followed a hierarchical measurement scheme to derive aboveground net primary productivity, functional attributes, and components of species diversity for grassland metacommunities. (a) Cover by species and aboveground biomass were measured in 76 cm diameter plots, (b) plots were randomly positioned in each 7 m diameter patch used for measurements to calculate patch-scale biomass, leaf dry matter content (LDMC), and diversity. (c) Patches were positioned in eight stands of restored grassland spanning two soil types (n = 52 stands per soil type). Soil types are demarcated by the line extending from upland at the left to lowland at the right of (c). Aboveground biomass and community LDMC per patch were calculated from measurements of canopy reflectance. (d) Adjacent patches were paired to create communities. Metacommunities were created from all possible pairs of communities. Alpha diversity of communities was calculated using variation in canopy reflectance between the two patches per community. Species dissimilarity between communities in each metacommunity was calculated using the between-community variation in reflectance. Illustrations of physical parameters are not scaled to actual size. ❖ www.esajournals.org 5 July 2020 ❖ Volume 11(7) ❖ Article e03178 communities that were located on each soil type (n = 24 communities per soil type). Metacommunities were created from all pairwise combinations of communities (n = 1128).
Metacommunity stability, diversity, and community attributes Species diversity and community attributes of communities were derived from remote measurements of canopy reflectance. Reflectance from each patch, including those for which cover was assessed, was measured at approximately 2week intervals during the 2016-2019 growing seasons (March-November). Diversity and community traits were calculated using reflectance data from the date each spring on which the across-patch mean of the enhanced vegetation index (EVI) peaked.
where NIR, Red, and Blue are reflectance in nearinfrared, red, and blue wavelengths, respectively, L is an adjustment for canopy background, G is a gain factor, and coefficients C1 and C2 correct for the influence of aerosols (L = 1, C1 = 6, C2 = 7.5, G = 2.5). Canopy reflectance was measured from an aerial platform (unmanned aerial vehicle [UAV], S1000; DJI, Shenzhen, China) using an ASD HandHeld2 spectroradiometer (spectral range of 350-1050 nm; ASD, Boulder, Colorado, USA). We measured reflectance by flying the GPSguided, rotary-wing UAV to a stationary position at 15.8 m height (25°field of view) above each patch. Reflectance was measured within 2 h of solar noon on cloudless days. All measurements were referenced to a Spectralon white reference panel at~15-min intervals. Reflectance was calculated by dividing radiance reflected from the plant canopy by radiance incident on the canopy. Incident radiation was considered the radiant flux reflected from a Spectralon white reference panel exposed to full sunlight. We used reflectance values from 5-nm intervals over the 660-760 nm range and 10-nm intervals over the remainder of the 350-1050 nm spectrum (80 wavebands) to calculate patch-scale diversity and community LDMC (see: Derivation of LDMC, α D , and dissimilarity from remote measurements). Reflectance for each patch and measurement date was normalized to a brightness of 1 (brightness normalization) to reduce the influence of canopy shading on model fits (Feilhauer et al. 2010).
The spring peak in aboveground biomass per patch (g/m 2 ; assumed equal to early-season ANPP) was calculated for each year as an exponential function of the spring maximum in EVI (adj. r 2 = 0.83, P < 0.0001; Appendix S1: Fig. S1). The biomass-EVI regression was developed using measurements of EVI and biomass from 32 patches. For each year, we summed calculated ANPP values of the two patches per community and two communities per metacommunity to derive ANPP for the community and metacommunity, respectively. Community LDMC per patch was estimated using a regression model developed from LDMC measurements on 52 single species stands (Polley et al. 2020).
We calculated metacommunity temporal stability (γ S ) and its spatial components, spatial asynchrony (ω) and alpha stability (α S ), following S. Wang, M. Loreau, C. de Mazancourt et al. (unpublished manuscript). Spatial asynchrony is defined as follows: where v ii is the sum of squares of community ANPP among years and v ij is the temporal covariance in ANPP between communities i and j (S. Wang, M. Loreau, C. de Mazancourt et al., unpublished manuscript). ω varies from 1 (perfect synchrony) to infinity (perfect asynchrony). The γ S is defined as the reciprocal of the squared coefficient of variation [1/CV 2 = 1/(σ/µ) 2 ] of metacommunity ANPP over 4 yr. Alpha stability is the square of the reciprocal of the weighted (by ANPP) average of the community CV in ANPP. As defined, γ S = α S × ω or, equivalently, log 10 γ S = log 10 α S + log 10 ω (S. Wang, M. Loreau, C. de Mazancourt et al., unpublished manuscript).
Derivation of LDMC, α D , and dissimilarity from remote measurements Partial least squares regression (PLSR) analysis was used to calculate LDMC per patch and two components of diversity, α D per community and species dissimilarity between communities, at the spring peak in EVI. Partial least squares ❖ www.esajournals.org 6 July 2020 ❖ Volume 11(7) ❖ Article e03178 regression identified a small set of uncorrelated latent variables using data from all wavebands. We used a cross-validation procedure to determine the number of latent variables required to predict the variable of interest (community LDMC, diversity). For each variable predicted, a PLSR model was fit to a reduced data set composed of every 7th observation of the variable and associated reflectance, beginning with the first entry in the full calibration data set, while minimizing the prediction error for unfitted data.
The process was repeated iteratively by beginning the split-sample procedure with data from each successive entry in the full calibration data set. For each variable predicted, we chose the PLSR model with the least number of latent variables for which the residuals from predictions of unfitted data (root mean predicted residual sum of squares; PRESS) did not differ significantly from the model with minimum PRESS. PLSR was fit using SAS 9.4 (SAS Institute, Cary, NC, USA). The LDMC per patch was estimated from a PLSR model developed previously (Polley et al. 2020). The PLSR model explained 73% of the variance in LDMC (%) in the calibration data. We averaged calculated LDMC values of the two patches per community to derive LDMC for each community.
We modeled α D of each 2-patch community as a function of the between-patch CV in reflectance at each waveband measured, as described by Polley et al. (2019b). The PLSR model of diversity was developed by using the CV in reflectance between patch pairs that were randomly selected (with replacement) from the total of 32 patches for which species diversity was measured (n = 48 patch pairs). The PLSR model explained 75% of the variance in α D of the calibration data set (Appendix S1: Fig. S2). The calibrated PLSR model was used to calculate α D for each community (n = 48) at the spring EVI peak each year (2016-2019).
Species dissimilarity between communities was calculated with a PLSR model developed using data from the 32 patches for which both patch-scale reflectance and cover by species were measured. We randomly paired patches (without replacement) to create 16 simulated communities, then calculated species dissimilarity between all pairwise groupings of communities (n = 120). This process was repeated to create a total of 240 pairwise comparisons of simulated communities. Species dissimilarity between communities was calculated using the abundancebased Bray-Curtis index. Species dissimilarity was modeled with PLSR using the ln-transformation of the between-community CV in brightness normalized values of reflectance at each of 80 wavebands (350-1050 nm range). This approach is based on the premise that between-community dissimilarity in species composition is related to between-community differences in reflectance. We used the calibrated PLSR model to calculate dissimilarity for all pairwise comparisons of communities measured in restored grassland each spring (n = 1128 pairwise comparisons).

Structural equation modeling
Structural equation models (SEMs) were used to evaluate causal relationships among variables that together affected ω between and α S of communities (Shipley 2000, Grace 2006). Models were fit using IBM SPSS AMOS 21 software. Asynchrony (log 10 ω) and alpha stability (log 10 α S ) additively account for 100% of variability in metacommunity stability (log 10 γ S ). Our primary objective in SEMs therefore was to determine contributions of diversity components and community functional traits to log 10 ω and log 10 α S , the components of log 10 γ S . The initial model of stability included five predictor variables. The model included two diversity metrics, α D and species dissimilarity in space, two community attributes, mean differences in biomass and community LDMC between communities, and one spatial variable, the Euclidean distance (m) between the mid-points of communities in each metacommunity (log e distance). Distance was included to account for spatial autocorrelation between communities. Primary predictors of log 10 ω were dissimilarity and community attributes. The α S was modeled initially as a function of α D . Communities that comprised each metacommunity were located on either the same soil type (silty clay or clay) or different soil types. We fit a separate model for each combination of soils (silty clay, clay, different soils). We iteratively removed parameters or pathways from the initial model in order to identify models with the most parsimonious fit. Parsimony was indicated by smaller values of the root mean square error of ❖ www.esajournals.org 7 July 2020 ❖ Volume 11(7) ❖ Article e03178 approximation (RMSEA) and Akaike's information criterion (AIC).

Diversity
Partial least squares regression accounted for 73% of the variance in species dissimilarity between communities created using data from patches on which species composition was measured (Fig. 4). Compositional differences between communities (dissimilarity) were associated largely with differences in reflectance in green (530-570 nm), red (620-700), far-red (725-780), and near-infrared (1030-1050) wavebands, as evidenced by departures in standardized weightings of PLSR coefficients from the zero line over these spectral ranges. Dissimilarity was positively correlated with α D (not shown; adj. r 2 = 0.06, P < 0.0001; α D range = 13.9-17.4) and with community differences in the attributes LDMC and biomass (not shown; adj. r 2 = 0.24 and 0.28, respectively, P < 0.0001).
Dissimilarity and α D explained little of the variance in temporal stability of grassland metacommunities. Dissimilarity explained 8% of variance in log 10 ω (Fig. 5). α D accounted for just 2% of variance in α S (P = 0.003; not shown). Dissimilarity-log 10 ω relationships were stronger when segregated by soil type. Dissimilarity explained 11% and 16% of variance in log 10 ω on silty clay and clay soils, respectively (not shown). α D accounted for 10% of variance in log 10 α S on the clay soil (not shown). Asynchrony was not correlated with either α D or community differences in α D across all 1128 pairwise comparisons of communities (P = 0.16 and 0.22, respectively).
Asynchrony was more highly correlated with community biomass and LDMC differences between communities than with species dissimilarly. Biomass and community LDMC differences accounted for 24% and 16%, respectively, of variance in log 10 ω of metacommunities from all soils combined (Fig. 5).

Structural equation models
Structural equation models were used to test interactive effects of diversity and community functional traits on the ω and α S of metacommunities. Metacommunities consisted of paired communities located on either the same soil type (silty clay or clay) or different soil types. One model adequately fit data for metacommunities Fig. 3. Bivariate relationships between temporal stability of aboveground net primary productivity (ANPP) of grassland metacommunities in spring and spatial asynchrony between (a) and average temporal stability of the ANPP of component communities (b; alpha stability). Lines represent linear regression fits to data (adj. r 2 = 0.28 and 0.58, respectively; P < 0.0001, n = 1128).
❖ www.esajournals.org 8 July 2020 ❖ Volume 11(7) ❖ Article e03178 composed of communities located on the clay soil combined with those located on different soil types (Fig. 6). A second, more complex model was required to fit data from the silty clay soil. Both SEMs confirmed the stronger correlation between γ S and α S than ω that was evident in bivariable relationships between metacommunity stability and its spatial components (Fig. 3). The α S in both models correlated more highly with differences in α D between communities than the weighted average of α D across communities comprising metacommunities.
Asynchrony and α S of metacommunities depended at least as much on biomass and LDMC differences between communities as on diversity differences alone (Fig. 6). Metrics of alpha and beta diversity (α D and dissimilarity) explained just 1.6% and 7.7% of the variance in Fig. 4. The relationship between species dissimilarity (Bray-Curtis index) measured between communities, each of two grassland patches, and species dissimilarity predicted with partial least squares regression (PLSR) using brightness normalized values of canopy reflectance (a; 5 latent variables; PRESS = 0.62; n = 240). The line is the 1:1 relationship with measured dissimilarity. Standardized weightings of regression coefficients for the PLSR model of species dissimilarity (b). Fig. 5. Bivariant relationships between spatial asynchrony in spring aboveground net primary productivity and mean values of species dissimilarity (a), biomass differences (b), and community leaf dry matter content (LDMC) differences (c) between grassland communities over 4 yr. Communities each consisted of two grassland patches. Solid lines are linear regression fits (n = 1128).
❖ www.esajournals.org 9 July 2020 ❖ Volume 11(7) ❖ Article e03178 α S and ω, respectively, in SEMs for metacommunities on the clay soil plus different soil types combined (χ 2 = 4.3; P = 0.89; df = 9; RMSEA < 0.0001, AIC = 28.3, n = 852). Substituting LDMC differences between communities for dissimilarity yielded a slightly less parsimonious model (RMSEA = 0.011, AIC = 34.0) that more than doubled the amount of variance in ω explained by the model (16.5%). Including both dissimilarity and LDMC differences between Fig. 6. Structural equation models of the roles of species diversity and community functional traits in regulating the spatial asynchrony and alpha stability of metacommunities. Metacommunities consisted of paired communities each of two 7 m diameter patches in restored grassland. Communities in each metacommunity were located on the same or different soil types (clay, silty clay). Community traits include differences in mean biomass and community leaf dry matter content (LDMC) between communities over 4 yr. The distance between communities in each metacommunity was included in models to account for spatial autocorrelation. Models were fit to data from communities located on the clay soil plus different soils combined or the silty clay soil only. Solid lines represent positive correlations. Dashed lines represent negative correlations. Standardized coefficients are listed beside the paths linking variables.
❖ www.esajournals.org 10 July 2020 ❖ Volume 11(7) ❖ Article e03178 units in the model only slightly increased the variance in ω explained (17.5%; χ 2 = 15.3; P = 0.22; df = 12; RMSEA = 0.018, AIC = 47.3). Leaf dry matter content differences were positively correlated with dissimilarity (r 2 = 0.22), but ω depended primarily on LDMC differences. Temporal stability on the silty clay soil also depended mostly on community functional attributes. A model that included dissimilarity with mediating effects on ω through biomass and LDMC differences between communities adequately fit data from the silty clay soil (χ 2 = 20.6; P = 0.11; df = 14; RMSEA = 0.04, AIC = 64.6, n = 276). Dissimilarity explained 29% and 49% of variance in LDMC and biomass differences, respectively. However, the most parsimonious model included biomass and LDMC differences alone. Combined, independent variables in this model for the silty clay soil explained 31% and 21% of the variance in ω and α S , respectively. Results from models for the soil groups considered confirmed that LDMC and biomass differences were positively correlated with dissimilarity. However, ω in models for both the silty clay soil and clay soil plus different soil types depended more on attribute differences than on diversity or diversity differences between communities.
Standardized regression coefficients linking diversity to stability and community traits to ω all were positive regardless of soil type. Asynchrony thus was increased by increasing both species dissimilarity and community differences in biomass and LDMC. Similarly, α S was increased by increasing community differences in α D .
Interestingly, biomass and LDMC differences between communities in metacommunities influenced α S as well as ω on the silty clay soil (Fig. 6). Standardized regression coefficients linking community traits to diversity were opposite in sign for α S and ω. Increasing biomass and LDMC differences increased γ S by increasing ω, but reduced γ S by reducing α D differences between communities and correlated α S . The net effect of differences in biomass and LDMC on γ S was negative for the silty clay soil. The standardized total effect of biomass differences on γ S was >7 times that of LDMC differences (−0.154 vs. −0.024).
Spatial autocorrelation between communities was much stronger for the silty clay than other soils, as evidenced by the greater influence of distance in the model for the silty clay soil (Fig. 6). All direct effects of distance (log e distance) were positive in both models, indicating that diversity components (dissimilarity, community differences in α D ) and biomass and LDMC differences between communities all were greater when communities comprising metacommunities were distant than near in proximity.

DISCUSSION
We analyzed effects of two components of species diversity (α D and species dissimilarity) together with spatial differences in two community functional traits (mean aboveground biomass and community LDMC) on temporal stability of spring ANPP of grassland metacommunities. Results demonstrated that productivity at the larger spatial scale of the metacommunity was stabilized more by spatial differences in community functional traits rather than by spatial differences in α D and greater dissimilarity between communities alone. Differences in community biomass and LDMC stabilized metacommunity productivity by increasing differences in the productivity responses of spatially distinct communities to interannual variability in precipitation (by increasing spatial asynchrony, ω), but de-stabilized ANPP on one soil type by reducing the temporal stability of communities (α S ). Results are consistent with the view that metacommunity stability depends on community differences in functional traits that couple plant growth to changes in resource availability (Polley and Wilsey 2018). Wilcox et al. (2017), for example, found that spatial synchrony was better correlated with synchrony in dynamics of spatially distinct populations within a species than with β diversity.
Diversity was estimated from measurements of canopy reflectance. Both α D and species dissimilarity between communities were well-predicted using spectral dissimilarities between patches in communities. Dalmayne et al. (2013) used Euclidean differences in reflectance per waveband to estimate species dissimilarity. We used the between-sample CV in reflectance per waveband to model species dissimilarity between communities, each composed of two patches. There are at least two explanations for ❖ www.esajournals.org 11 July 2020 ❖ Volume 11(7) ❖ Article e03178 the observed correlation between spectral and species dissimilarities. First, the link between spatial differences in canopy optical properties and species composition may be indirect. The spatial differences in canopy structure or chemistry that were detected by remote sensing may reflect spatial differences in soils or disturbance that are in turn associated with differences in species composition (Rocchini et al. 2010, Dalmayne et al. 2013, Möckel et al. 2016. Secondly, spatial differences in canopy optical properties may reflect structural or chemical differences among species themselves (Polley et al. 2019b). This latter explanation for an association between spectral and species dissimilarities is the more plausible in this study given that the grassland sampled was both rather small and uniformly managed. More generally, results provide evidence that remote sensing is useful for assessing spatial components of species diversity. Grassland temporal stability correlated with differences in community attributes that led to differences in biomass production in space and operated mainly through spatial asynchrony. Asynchrony in growth has been shown to be greater among functionally dissimilar species at local scales (Roscher et al. 2011, Polley et al. 2013, Craven et al. 2018, Lepš et al. 2018). Our results demonstrate that asynchrony is linked to differences in functional properties between communities.
Communities differed in attribute values partly because they differed in composition (exhibited heighted species dissimilarity), but dissimilarity explained <50% of the variation in community differences in mean biomass and LDMC. One implication is that selection for spatial differences in dominant trait values and the expression of traits in biomass production proceeded partly independently of taxonomic variation. This result is perhaps not surprising, as species dissimilarity and community attributes reflect different aspects of biodiversity. Dissimilarity may increase temporal stability by favoring spatial differences in properties such as phenology and rooting depth (Fargione and Tilman 2005, Dornbush and Wilsey 2010) that are not closely associated with community-weighted values of leaf or other growth-related traits. Conversely, community attributes are an index of functional composition. Communities that differ in functional composition may differ in biomass resistance to and recovery following precipitation change and other fluctuations.
Attribute effects on temporal stability depended on soil type. Community differences in LDMC and biomass increased γ S on the clay soil and clay plus silty clay (different) soils combined, but decreased γ S on the silty clay soil. The latter resulted because α S was negatively correlated with community differences in functional attributes. This divergence in functional regulation between soils implies that attribute effects on stability were mediated by biotic or abiotic effects associated with soil type. ANPP is highly responsive to precipitation in this grassland (Polley et al. 2020). ANPP, thus, likely is sensitive to soil effects on water availability. Water-holding capacity is lower on the silty clay than clay soil (Polley et al. 2019a). It is plausible therefore that communities on the silty clay soil experienced greater or more frequent water limitation which led, in turn, to a stronger impact of growth-related attributes on temporal stability on this soil type. S. Wang, M. Loreau, C. de Mazancourt et al. (unpublished manuscript) also reported an abiotic dependence of diversity-stability relationships. Precipitation and other climatic variables influenced the relationship between β D and spatial asynchrony of simulated metacommunities in their study.
Metacommunities were stabilized more by α S than by ω. A similar trend has been reported in other systems (Wilcox et al. 2017;S. Wang, M. Loreau, C. de Mazancourt et al., unpublished manuscript). The differences in standardized effects of stability components on γ S (0.78-0.87 and 0.61-0.64 for α S and ω, respectively) as derived from SEMs were, nevertheless, relatively small for a grassland that had been seeded with a defined species mixture 6 yr prior to our measurements. The significant ω effect on temporal stability reflects large spatial differences in ANPP responses to precipitation. Temporal variation in ANPP, in turn, reflects spatial variation in the expression of functional attributes, which likely resulted partly from spatial variation in establishment and growth of both seeded perennials and a cohort of invasive annual species. Spatial variation in establishment likely increased species dissimilarity in space. The imputed importance of spatial variation in species establishment ❖ www.esajournals.org and sorting is consistent with the contribution of spatial autocorrelation (distance) to diversity and temporal stability in SEMs.
The components of biodiversity that we studied explained relatively little of the variance (<32%) in spatial components of temporal stability, particularly for metacommunities on clay plus different soil types combined. Prediction may have been improved had additional components of biodiversity been measured and included in the model (Craven et al. 2018). Inclusion of functional diversity (diversity calculated using species traits), for example, may have improved model explanation of stability components, especially if derived using multiple species traits. Unfortunately, few data exist for some of the traits potentially associated with species niche differences, including phenology and rooting depth.
Our results indicate that diversity had limited direct influence on temporal stability of restored grassland. Metacommunity stability was more directly correlated with community differences in functional attributes that link plant growth to environmental fluctuations. Results imply that temporal stability of ANPP can be improved in some grasslands by practices that increase spatial differences in the dominance of plant attributes that couple growth and resource availability. We recommend increased attention to developing techniques to monitor species-abundance weighted functional traits that influence productivity and its response to the environment.