Stability and asynchrony of local communities but less so diversity increase regional stability of Inner Mongolian grassland

Extending knowledge on ecosystem stability to larger spatial scales is urgently needed because present local-scale studies are generally ineffective in guiding management and conservation decisions of an entire region with diverse plant communities. We investigated stability of plant productivity across spatial scales and hierarchical levels of organization and analyzed impacts of dominant species, species diversity, and climatic factors using a multisite survey of Inner Mongolian grassland. We found that regional stability across distant local communities was related to stability and asynchrony of local communities. Using only dominant instead of all-species dynamics explained regional stability almost equally well. The diversity of all or only dominant species had comparatively weak effects on stability and synchrony, whereas a lower mean and higher variation of precipitation destabilized regional and local communities by reducing population stability and synchronizing species dynamics. We demonstrate that, for semi-arid temperate grassland with highly uneven species abundances, the stability of regional communities is increased by stability and asynchrony of local communities and these are more affected by climate rather than species diversity. Reduced amounts and increased variation of precipitation in the future may compromise the sustainable provision of ecosystem services to human well-being in this region.


Introduction
The ability of ecosystems to stably provide biological products and services such as biomass production for human well-being Tilman et al., 2014, Tilman et al., 2006 is being threatened by species loss (Cardinale et al., 2012;Harrison et al., 2015;Isbell et al., 2017; Box 1. Theory and glossary. Theoretically, the local population coefficient of variation (CV, i.e., the ratio of mean to standard deviation or the inverse of stability) can be stepwise upscaled to the regional community CV via either the local community CV (pathway I) or the regional population CV (pathway II) (Box 1-figure 1). In each step, synchrony (inverse of asynchrony) measures the proportion of CV upscaled from the lower to the higher hierarchical level. Specifically, synchrony takes value between 0 (perfectly asynchronous) and 1 (perfectly synchronous) and the CV at the higher level is the product of synchrony and the CV at the lower level (Loreau and de Mazancourt, 2008;Thibaut and Connolly, 2013;Wang et al., 2019). Along pathway I, the local population CV first upscales with local species synchrony to the local community CV (step A); then, the local community CV upscales with regional community synchrony to the regional community CV (step B). Along pathway II, the local population CV first upscales with regional population synchrony to the regional population CV (step A); then the regional population CV upscales with regional species synchrony to the regional community CV (step B). In this study, we use superscripts P and C to designate hierarchical components at population and community levels, respectively, and superscripts L and R to designate local and regional scales, respectively. We use superscript P→C and L→R to designate scaling up from populations to communities and from local to regional measures, respectively. All measures were estimated with all species or only dominant species, the latter designated with subscript d.
Temporal CV and synchrony Local population CV (CV P,L ): Defined as the weighted-average local population CV of plant aboveground biomass across species and local communities. Hypothesis: positively or negatively influenced by alpha diversity (Thibaut and Connolly, 2013;Tilman et al., 2006). Local species synchrony (φ P→C,L ): Defined as the weighted-average synchronous biomass dynamics among local populations within local communities. Hypothesis: negatively influenced by alpha diversity (Loreau and de Mazancourt, 2008;Thibaut and Connolly, 2013). Local community CV (CV C,L = CV P,L × φ P→C,L ): Defined as the weighted-average local community CV of biomass among local communities. Hypothesis: negatively influenced by alpha diversity (Loreau and de Mazancourt, 2008;Thibaut and Connolly, 2013;Tilman et al., 2014). Regional community synchrony (φ C,L→R ): Defined as the weighted-average spatial synchronous biomass dynamics among local communities. Hypothesis: negatively influenced by beta diversity .
Regional population synchrony (φ P,L→R ): Defined as the weighted-average spatial synchronous biomass dynamics among local populations of the same species. Hypothesis: negatively influenced by spatial heterogeneity . Regional population CV (CV P,R = CV P,L × φ P,L→R ): Defined as the weighted-average regional population CV of biomass across species. Hypothesis: positively or negatively influenced by gamma diversity . Regional species synchrony (φ P→C,R ): Synchronous biomass dynamics among regional populations of different species. Hypothesis: negatively influenced by gamma diversity . Regional community CV (CV C,R ): Defined as the CV of regional community biomass. Can be upscaled via aggregating local communities (pathway I, CV C,R = CV C,L × φ C,L→R ) or organizing regional populations (pathway II, CV C,R = CV P,R × φ P→C,R ). Hypothesis: negatively influenced by gamma diversity .

Species diversity
Alpha species diversity: Local species richness (S α ) or effective species richness (D α ). Beta species diversity: Cross-locality dissimilarity of species richness (S β ) or effective species richness (D β ). Gamma species diversity: Regional species richness (S γ ) or effective species richness (D γ ).

Climatic variables
Mean annual precipitation (MAP): Cross-site averaged mean annual precipitation (or annual temperature, MAT). Local and regional precipitation CVs (CV P L and CV P R ): Measuring precipitation (or temperature, CV T L and CV T R ) variability with its local and regional CVs. Regional precipitation synchrony (φ P L→R ): Spatial synchronous dynamics of precipitation (or temperature, φ T L→R ).
in China across an area of >166,894 km 2 and monitored the yearly dynamics of productivity over five consecutive years (Figure 1a). The Inner Mongolian grassland represents a typical part of the Eurasian grassland biome and is crucial in providing biological products and services to human societies living there Kang et al., 2007). In this region, plant community productivity and species richness and composition are driven by climatic factors, that is, temperature and precipitation (Bai et al., 2004;Hu et al., 2018;Ma et al., 2010;Wang et al., 2020;Xu et al., 2015). These have changed considerably during the past decades (Huang et al., 2015;Piao et al., 2010) with largely unknown ecological consequences, especially at large spatial scales.
To facilitate the large-scale stability investigation, we employed a simulated landscape method (Hautier et al., 2018;van der Plas et al., 2016) to construct large-scale, that is, regional communities consisting of two local communities (two observed sites) separated by 17-987 km ( Figure 1a). Briefly, each regional community was constructed by randomly choosing two distant local communities without replacement (to ensure replicate regional communities were not sharing local communities; Figure 1b). We did not consider scenarios including more than two local communities in each regional community because the resulting small number of replicate regional communities would have prevented a statistical analysis (but see Appendix 1-figure 3 for a three-local-community scenario). Based on the above framework, we investigated how asynchronous population dynamics among species, in particular dominant ones, at local and regional scale contributed to the stability of local and regional communities in the study region. We also tested how local and regional community dynamics were driven by species diversity or affected by climatic factors such as precipitation and its temporal variation. First, we analyzed stability variables with general linear models (GLMs) to identify important relationships. We then used this information together with theoretical considerations to construct path-analytic diagrams from structural equation models (SEMs) relating the regional community coefficient of variation (CVs; i.e., inverse of stability) of plant aboveground biomass to its hierarchical (a) (b) Random sample Random sample Regional community (randomly sampled) Figure 1. Geographic distribution of surveyed sites with site numbers (a) and a simplified case (seven-site) for illustrating construction of regional communities aggregating two local communities (b). In (a), red circles represent sites included in constructing regional communities. (b) shows a simplified case illustrating the construction of two sets of regional communities using random sampling without replacement to ensure the regional communities within each set do not share local communities.
The online version of this article includes the following figure supplement(s) for figure 1: Figure supplement 1. Regional synchronies of temperature (a) and precipitation (b), all-species measure of regional community coefficient of variation (CV, inverse of stability, c) and all-species (d-f) and dominant-species (g-i) measures of regional community synchrony (inverse of asynchrony), regional population CV, and regional population synchrony in relation to distance.
components and all-or dominant-species diversity as well as climatic factors (Thibaut and Connolly, 2013;Tilman et al., 2014;Wang et al., 2019;Wang and Loreau, 2016;Wang and Loreau, 2014). As commonly done in such studies, all analyses were conducted using inverse values of stability and asynchrony, that is, CV and synchrony. We present an overview of the upscaling models and a glossary of terms in Box 1.

Results
Part I: Analysis using all species We first analyzed variation in the regional community CV in relation to its hierarchical components including all species and in relation to all-species diversity indices as well as climatic factors. We found that the regional community CV was positively associated with the local community CV and regional community synchrony (step B of upscaling pathway I; Figure 2a and b, Figure 3a). The local community CV in turn was positively related to the local population CV and local species synchrony (step A in upscaling pathways I; Figure 2e and f, Figure 3a). Along the upscaling pathway I, the CVs decreased from 0.76 for the local population CV to 0.38 for the local community CV and further to 0.29 for the regional community CV (Figure 3a), as a result of a lower local species synchrony (mean = 0.49) compared with regional community synchrony (mean = 0.78; Figure 3a). Alternatively, the regional community CV was positively related to the regional population CV and regional species synchrony (step B of upscaling pathway II; Figure 2c and d, Figure 3a). The regional population CV was in turn positively related to the local population CV but not to regional population synchrony (step B of upscaling pathway II; Figure 2g and h). However, the path from regional population synchrony to regional population CV was suggested by theory and therefore included in the SEM, where it became significant (Figure 3a). Along the upscaling pathway II, the CVs declined from 0.76 for the local population CV to 0.71 for the regional population CV and further to 0.29 for the regional community CV (Figure 3a), as a result of a higher regional population synchrony (mean = 0.94) compared with regional species synchrony (mean = 0.41; Figure 3a). Figure 2. The regional community (a-d), local community (e-f), and regional population (g-h) coefficients of variation (CVs) in relation to their hierarchical components using all-species measures. Solid black lines represent significant (p<0.05) and marginally significant (p<0.10) relationships, and dashed gray lines represent nonsignificant (p>0.10) relationships (see 'Materials and methods' for details and Box 1 for glossary). Thin grey lines represent Figure supplement 1. All-species estimate of regional community coefficient of variation (CV, inverse of stability) in relation to dominant-species estimates of local community CV (a), regional community synchrony (inverse of asynchrony, b), regional population CV (c), and synchrony (d) as well as dominant-species estimates of local community CV (e, f) and regional population CV (g, h) in relation to their hierarchical components.

Figure 2 continued
We found that all-species diversity indices had relatively weak impacts on CVs and synchronies across ecological organization levels (Figure 4a and c-e; see 'Materials and methods' for calculating species diversity indices across scales). Although correlation (Appendix 1-figure 2a) and regression ( Figure 4b) analyses showed that the regional population CV was positively related to gamma diversity, this was not supported by the final path analysis model ( Figure 3a). However, the local population CV was positively and the local species synchrony negatively related to alpha diversity (Figure 3a, Figure 4f and g).
Local species synchrony increased with the local precipitation CV but no relations between the regional community CV or its other components and climatic variation were found (Figure 3a, Figure 4h and i).

Part II: Analysis using only dominant species
Considering only dominant species in the hierarchical components was sufficient to explain a large amount of variation in regional community CV. For the upscaling pathway I, the regional community CV was positively related to local community CV and regional community synchrony (step B), with the explanatory power reduced to 0.71 from 0.98 for the analysis using all species (comparison of Figure 3b with Figure 3a). The local community CV in turn was positively related to the local population CV and to local species synchrony (step A; Figure 3b). For the upscaling pathway II, the regional community CV was positively related to the regional population CV and regional species synchrony (step B), with the explanatory power reduced to 0.69 from 0.98 for the analysis using all species. The regional population CV was positively related to the local population CV and regional population synchrony. Dominant species as a group explained more than half of the variance of CVs and synchronies estimated with all species (explanatory power, R2 ,>0.50, Figure 3-figure supplement 1), except for the regional population synchrony ( R2 = 0.14, Figure 3-figure supplement 1h).
Similar to the analysis with all species, diversity indices of only dominant-species had relatively weak impacts on CVs and synchronies across organizational levels (Figure 4-figure supplement 1a-g). Although correlation analyses showed that the regional population CV was positively related to gamma diversity (Appendix 1-figure 2b), this was again not supported by the final path analysis model ( Figure 3b). However, for dominant species, the local population CV was increased by alpha diversity and reduced by larger mean values of precipitation (

Discussion
Based on a multiyear region-scale survey in Inner Mongolian grassland, we investigated stability (inverse of CV) and asynchrony (inverse of synchrony) across spatial scales and analyzed influences of species diversity, dominant species, and climatic factors on them. We found that the regional community stability was related to the stability and asynchronous dynamics of local communities. In addition, stability and asynchrony were -albeit weakly -impacted by species diversity. Compared with the dynamics of all species, the dynamics of only dominant species had also good predictive power, indicating that these species were important drivers of grassland stability in the region. Furthermore, decreasing mean and increasing interannual fluctuation of precipitation could, respectively, destabilize dominant species and synchronize population dynamics within local communities, impairing stability at the regional scale. . Path analysis models relating the regional community coefficient of variation (CV, all-species measure) to its hierarchical components and species diversity indices estimated with all species (a, the mean values of CVs and synchronies are in brackets) and only dominant species (b) as well as climatic factors. These diagrams combine different upscaling pathways (pathway I, left side with red arrows; pathway II, right side with blue arrows). Solid and dashed arrows, respectively, represent positive and negative paths, and numbers near arrows are standardized path coefficients. The significance level of each path is indicated by * when p<0.05 or # when p<0.10 (see 'Materials and methods' for details and Box 1 for glossary). See Figure supplement 1. All-species estimates (vertical axes) of coefficients of variation (CVs, inverse of stability) and synchronies (inverse of asynchrony) across hierarchical levels of ecological organization in relation to their dominant-species counterparts (horizontal axes) (a-b, regional community CV; c, regional community synchrony; d, regional species synchrony; e, local community CV; f, regional population CV; g, local species synchrony; h, regional population synchrony; i, local population CV). 4. Regional community coefficient of variation (CV) (a), regional population CV (b), regional species synchrony (c), regional community synchrony (d), local community CV (e), local population CV (f), and local species synchrony (g) in relation to species diversity (effective species richness) as well as local species synchrony and local population CV, respectively, in relation to local precipitation variability (h) and precipitation (i) using all-species measures. Solid black lines represent significant (p<0.05) and marginally significant (p<0.10) relationships, and dashed gray lines represent nonsignificant (p>0.10) relationships (see 'Materials and methods' for details and Box 1 for glossary). Thin grey lines represent relationships for 1000 sets of resampled regional communities (n=10 for each set). See Figure supplement 1. Regional community coefficient of variation (CV, inverse of stability) estimated with all species (a) and dominant-species estimates of regional population CV (b), regional species synchrony (c), regional community synchrony (d), local community CV (e), local population CV (f), and local species synchrony (g) in relation to species diversity (effective species richness) as well as dominant-species estimates of local species synchrony and local population CV in relation to local precipitation variability (h) and precipitation (i), respectively.

Stability across ecological hierarchies
We investigated ecosystem stability across ecological hierarchies with two alternative upscaling pathways , and both of them showed gradually increasing stability from low to high organization levels due to species insurance effects and spatial insurance effects of populations and communities, caused by asynchronous dynamics among species and localities ( Figure 3a). These patterns are consistent with recent studies carried out at single sites constructing multiple adjacent plots within meta-communities in grassland ecosystems (Hautier et al., 2020;McGranahan et al., 2016;Wang et al., 2019;Wilcox et al., 2017;Zhang et al., 2019) and at the regional scale in marine ecosystems Thorson et al., 2018), as well as recent theoretically proposed positive invariability-area relationships Wang et al., 2017). These results suggest that, at large spatial scales, spatial heterogeneity is important in maintaining stability; losing this heterogeneity (Fahrig et al., 2011;Gámez-Virués et al., 2015) can impair stability.
We found that species insurance effects (local species asynchrony; Tilman et al., 2014;Yachi and Loreau, 1999) were stronger in maintaining stability at the regional scale than the spatial insurance effects of distant populations and communities despite the large extent and thus expected spatial heterogeneity of our study region. This result is consistent with a recent investigation in marine plant communities  but different from that in fish communities (Thorson et al., 2018). In our study, the region-wide synchronous variation in precipitation (mean = 0.86, ranged from 0.62 to 1.00) (Figure 1-figure supplement 1b) potentially decreased the spatial heterogeneity and increased the relative importance of the species insurance effect. The regulation of spatial insurance effects on the stability of fish communities at regional scale may result from their high mobility. Fish can move toward their preferred environmental conditions, causing more variable spatial population patterns than those found in plants, potentially strengthening the spatial insurance effects. In plant communities, the strong species insurance effect suggests that regional community stability to a large part reflects the stability of local communities, which have predominantly been considered in previous studies (Ma et al., 2017;Tilman et al., 2006;Xu et al., 2015;Yang et al., 2012).

Influence of species diversity, dominant species, and precipitation on ecosystem stability
We only detected stabilizing impacts of species diversity at local but not at regional scale (Figure 3, Figure 4f and g). The observed negative species richness-local population stability relationship is in line with theoretical and experimental studies (Lehman and Tilman, 2000;Tilman, 1999;Tilman et al., 2014, Tilman et al., 2006, proposing that competition between coexisting species for resources in species-rich communities leads to low population stability. The observed positive species richnesslocal species asynchrony relationship is also expected by theory based on the higher probability of species-rich communities to contain species that are different in responding to environmental fluctuations (Tilman et al., 2014;Yachi and Loreau, 1999).
Previous studies reported mixed impacts of species diversity on stability and asynchrony at scales beyond local. Some studies found significant influences (e.g., Hautier et al., 2020;Liang et al., 2021;Qiao et al., 2022;Wang et al., 2019) and others found none (e.g., Wilcox et al., 2017;Yang et al., 2022;Zhang et al., 2019). It has been argued (Hautier et al., 2020) that investigations within a single site (Zhang et al., 2019) or multiple sites with nonstandardized experimental protocols (Wilcox et al., 2017) may mask stabilizing effects of species diversity at regional scale. However, even though this study used a multisite dataset across a large region with a standardized survey protocol, it still could not detect significant effects of species diversity at the regional scale ( Figure 3). The highly uneven distribution of species abundances could in part have been responsible for this . Nevertheless, even when only dominant species were considered, we could still not find the expected relationship between species diversity and regional asynchrony and stability.
It is conceivable that in natural grassland ecosystems, which are often characterized by high unevenness (Dee et al., 2019;Jiang et al., 2009;Smith and Knapp, 2003), species have co-evolved over time in such a way as to maintain high stability at a level of species diversity established over a longer time span. This observation has been made in a long-term grassland biodiversity experiment, where a history of co-occurrence led to greater community stability in response to a flooding event (van Moorsel et al., 2021) and repeated exposure to drought led to increased species complementarity in response to drought (Chen et al., 2022). The absence of strong species diversity-stability relationship at regional scales in observational studies is not the same as rapid biodiversity loss from an established level typically simulated in experiments. By extension, regional species loss in the study area due to, for example, land-use and climate change in the future may still pose a threat to regional grassland stability.
The strong influence of precipitation on the productivity of different species  may have weakened insurance effects of species diversity in this study. Strong fluctuations in precipitation in dry grassland may force similar responses among different species, decreasing the dissimilarity and thus compensatory dynamics among species. This speculation is supported by the low local species asynchrony under high precipitation fluctuation (Figure 3a, Figure 4h). Furthermore, for the dominant-species measures, we also found decreased local population stability under low precipitation means (Figure 3b, Figure 4-figure supplement 1i), potentially due to the decreasing mean-to-standard deviation ratio caused by the dominant-species biomass production being more steeply related to precipitation amount than to its standard deviation . The study region has been experiencing a pronounced decrease in precipitation and an increase in its variability during the past decades (Huang et al., 2015;Piao et al., 2010;Tao et al., 2015). Our results indicate that these changes in precipitation regimes may present a key threat to the sustainable provision of biological products and services to human well-being in the region.

Study region and plant community survey
The Inner Mongolian temperate grassland has a continental monsoon climate with a short and cool growing season (from May to October, averaged temperature 12.9-18.4°C across sites during the studied period from 2012 to 2016), concentrating ~90% of the annual precipitation (averaged precipitation 186.2-398.0 mm across sites from 2012 to 2016) . We established a 5-year (2012-2016) region-scale survey over this area (latitudes 39.34-49.96°N, longitudes 107.56-120.12°E), covering different grassland types (Figure 1a; Wang et al., 2020). There were 21 sites with 4-5 consecutive years of data. The sample plots of each site were randomly selected, excluding anthropogenic disturbances (e.g., overgrazing and heavy mowing). At each site, we marked three 1 × 1 m quadrats along the diagonal of a 10 × 10 m plot, harvested all living plant tissues, sorted them to species, and then oven-dried and weighed the harvested material to obtain aboveground biomass and species richness (for details, see Wang et al., 2020).

Construction of regional communities
We constructed regional communities consisting of two local communities with a simulated landscape method (Hautier et al., 2018;van der Plas et al., 2016). Specifically, the 21 local communities (sites) were randomly separated into 10 regional communities without replacement (two local communities for each regional community with one remainder) to ensure that they were independent between each other (see Figure 1b for a simplified seven-site case). We repeated this random resampling process 1000 times, resulting in 1000 resampled sets, each containing 10 regional communities that were independent of each other.
Temporal CV, synchrony, and species diversity across ecological hierarchies A regional community includes M local communities (M = 2 for the current case) and S species (see Box 1 and Appendix 1-table 1). Its temporal dynamics can be described with a matrix with elements u P,L (i, k) for the mean abundance of species k in locality i, and a matrix with elements v P,L (ij, kl) = cov(N-P,L (i, k, t), N P,L (j, l, t)) for the covariance between abundances of species k in locality i and species l in locality j over time t. Here, N denotes population abundance, and the superscripts indicate 'population' (P) and 'local' (L). To estimate CV and synchrony with only dominant species (designated with subscript d, see Appendix 1 for detailed mathematical derivation), we introduce a matrix d P with elements d P (I, k) set to 1 if the species k is a dominant species, otherwise, 0. We defined dominant species as those whose biomass contributed to >5% of the total biomass of the regional community  over the years of the survey (see Appendix 1-figure 1).
The local population CV can be upscaled to the regional community CV via two alterative pathways I or II (Box 1; Wang et al., 2019). First, we estimated the local population CV of all species (CV P,L , Equation 1a) and only dominant species (CV d P,L , Equation 1b) as follows: Here, u C,R is the average total biomass of the regional community over time and the superscripts indicate 'community' (C) and 'regional' (R). Second, we estimated the local species synchrony of all species ( Here, the superscript (P→C, L) indicates upscaling along pathway I, step A, from local population CV and local species synchrony to local community CV. Third, the local community CV of all species (CV C,L , Equation 3a) and only dominant species (CV d C,L , Equation 3b) was estimated: The regional community synchrony of all species (φ C,L→R , Equation 4a) and only dominant species (φ d C,L→R , Equation 4b) can be estimated as follows: The superscript (C, L→R) indicates upscaling from local community CV and regional community synchrony to regional community CV (pathway I, step B). Here, v C,L (ii) = ∑ kl v P,L ( ii, kl ) is the variance of community biomass at locality i over time. Finally, the regional community CV of all species (CV C,R , Equation 5a) and only dominant species (CV d_C C,R , Equation 5b) along the upscaling pathway was estimated: Note that the regional community CV of only dominant species is presented here only for completeness, but in our empirical analysis we were more interested in the relationship between the regional community CV of all species and dominant species dynamics. That is, we wanted to test the explanatory power of dominant-relative to all-species components in predicting the all-species regional community CV (i.e., CV C,R ; Box 1).
Along pathway II, local population CV and regional population synchrony scale up to regional population CV at first (pathway II, step A). The regional population synchrony of all species (φ P,L→R , Equation 6a) and only dominant species (φ d P,L→R , Equation 6b) is estimated as follows: The regional population CV of all species (CV P,R , Equation 7a) and only dominant species (CV d The regional population CV and regional species synchrony scale up to the regional community CV (pathway II, step B). The regional species synchrony of all species (φ P→C,R , Equation 8a) and only dominant species (φ d P→C,R , Equation 8b) is estimated as follows: Here, v P,R (kk) is the variance of population biomass of species k at the regional scale over time. The regional community CV of all species (CV C,R , Equation 9a) and only dominant species (CV d_P C,R , Equation 9b) according to the upscaling pathway II can be estimated as follows: Here again the regional community CV of only dominant species is listed for completeness, but the regional community CV of all species was related to dominant species dynamics to test the explanatory power of dominant-relative to all-species components in predicting the all-species regional community CV along upscaling pathway II (see Box 1).
The regional community CV estimated with the two alterative upscaling pathways is the same when using all-species measures, but can be slightly different when using dominant-species measures (see Appendix 1-1.6 for details), which is why we used two abbreviations, that is, CV d_C C,R and CV d_P C,R . We estimated two alternative species diversity indices across ecological hierarchies, species richness (S) and effective species richness (D). The alpha (S α ) and gamma species richness (S γ ) were defined as the total number of species at local and regional scales, respectively, and the beta species richness (S β =S γ / S α ) was used to measure dissimilarity among localities. Specifically, the alpha (S α ) and gamma (S γ ) species richness were estimated as multiple-year mean (S α ) and multiple-year pooled species number (S γ ) of the two local communities. To account for uneven species abundances in the study region, we also used effective species richness, the antilog of Shannon-Wiener diversity (D = e H' ), reflecting how many species with an even abundance distribution would produce the same Shannon-Wiener diversity as observed for the actual uneven community . The alpha (D α ) and gamma (D γ ) effective species richness thus represented the Shannon-Wiener diversity at local and regional scales, respectively, with beta effective species richness (D β = D γ /D α ) measuring its cross-locality dissimilarity. These species diversity indices were estimated with either all species or only dominant species.

Climatic data
Based on monthly climatic data collected from 119 climate stations and 2 km resolution digital elevation over this region, we calculated site-specific mean temperature and precipitation using the simple kriging method and spherical model of geostatistical analysis in ArcGIS software (Environmental Systems Research Institute Inc, Redlands, CA). We calculated mean annual temperature (MAT) and annual precipitation (MAP) (only monthly data from May to October were used as plants are active only in this period), as well as their CVs across spatial scales from 2012 to 2016. CVs of temperature and precipitation at the local (CV T L and CV P L ) and regional scales (CV T R and CV P R ), as well as their among-site synchronies (φ T L→R and φ P L→R ), were estimated with the methods used for local and regional community CVs and regional community synchrony.

Statistical analysis
We analyzed the influence of distance between spatially separated local communities (i.e., sites) within regional communities on regional synchronies of temperature and precipitation, regional community CV, as well as all-species and dominant-species measures of regional population CV and synchrony and regional community synchrony with linear regressions. However, spatial distance only influenced regional synchronies of temperature and precipitation (Figure 1-figure supplement 1), and thus distance was not considered further in subsequent analyses.
We used correlation analyses, linear regression analyses, and path analyses to investigate the regional community CV in relation to its hierarchical components, species diversity indices, and climatic factors. For the path analyses, we established models considering different upscaling pathways and different species diversity indices to explain variation in the all-species regional community CV either using all species or only dominant species. We started with initial models as close as possible to paths proposed in theoretical studies Wang and Loreau, 2016;Wang and Loreau, 2014). GLMs were used to analyze the proposed paths. Then, we modified the initial path-analysis models (Appendix 1- figure 4 and Appendix 1-figure 5) to eliminate nonsignificant paths until only significant or marginally significant paths remained or a minimal value of Akaike's information criterion for small sample size (AICc) was reached (Brewer et al., 2016). Path coefficients of the final models were quantified using the piecewiseSEM package (Lefcheck, 2016) of R 3.6.3 (R Development Core Team, 2013. We used a randomized examination method to investigate the statistical significance of the above analyses (Edgington and Onghena, 2007;Efron and Tibshirani, 1994). Specifically, considering the 10 independent regional communities per sampled set, all above statistical analyses were conducted within each set, resulting in 1000 statistics. Taking the correlation analysis as an example, we calculated the mean correlation coefficient ( − ρ ) of the 1000 sets and considered it to be statistically significant or marginally significant if the proportion of ρ < 0 (P -ρ ) (or ρ > 0, P +ρ ) was lower than 0.05 or 0.10 when − ρ > 0 (or − ρ < 0), respectively. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Appendix 1

Mathematical derivation for partitioning coefficients of variation and synchronies across ecological hierarchies into dominant and subdominant species groups
In the following, we introduce methods partitioning (temporal) coefficients of variation (CVs, inverse of (temporal) stability) and synchronies (inverse of asynchrony) across ecological hierarchies into dominant (relative species abundance >5%, see Appendix 1- figure 1 for details) and subdominant species groups. We do not repeat theoretical derivations relating CVs and synchronies across different hierarchies here because these can be found elsewhere (Loreau and de Mazancourt, 2008;Thibaut and Connolly, 2013;Wang et al., 2019;Wang and Loreau, 2016;Wang and Loreau, 2014). Symbols and descriptions used in the following partitions can be found in Box 1 and Appendix 1-table 1.
Appendix 1-table 1. Notation summary for climatic factors, species diversity indices, coefficients of variation (CVs, inverse of stability), and synchronies (inverse of asynchrony) across spatial scales and hierarchical levels of ecological organization.

Climatic factors
MAT Cross-site averaged mean annual temperature

MAP
Cross-site averaged mean annual precipitation Local CV of precipitation φ T

L→R
Regional temperature synchrony φ P

L→R
Regional precipitation synchrony CV T R Regional CV of temperature CV P ) or organizing regional populations (pathway II:

Regional CV of precipitation
We consider a regional community reached a stationary state, which includes M localities (e.g., sites or local communities) and S species. This regional community can be described with a matrix of (temporal) mean species abundance with elements u P,L (i, k), that is, the mean abundance of species k in locality i, and a (temporal) variance-covariance matrix of species abundances with elements v P,L (ij, kl) = cov (N P,L (i, k, t), N P,L (j, l, t)), that is, the covariance between abundances of species k in locality i and species l in locality j over time t. In addition, we introduce two matrixes, d P and s P , to represent the dominant and subdominant species of the regional community, respectively. For the d P , it has M rows and S columns, representing numbers of locality and species of the regional community, and has elements d P (i, k), that is, the kth species of the ith locality, which is set to 1 if the kth species is a dominant species at the regional scale, otherwise, 0. Similar procedure is used to conduct the s P , in which subdominant species are set to 1, otherwise, 0.

Partitioning local population CV into dominant and subdominant species groups
The local population CV (CV P,L ) is defined as the weighted-average local population CV (Thibaut and Connolly, 2013;Wang et al., 2019;Wang and Loreau, 2016;Wang and Loreau, 2014): We rewrite this equation with introduced two matrixes (d P (i, k) and s P (i, k)) to separate the local population CV (CV P,L ) into its dominant (CV d P,L ) and subdominant (CV s P,L ) species group components:

Partitioning local species synchrony into dominant and subdominant species groups
The local species synchrony (φ P→C,L ) is defined as the weighted-average synchronous dynamics among populations of different species within local communities Wang and Loreau, 2016;Wang and Loreau, 2014): where ω P→C,L (i) and φ P→C,L (i) are the contribution of local population variance of the ith community to the sum of variance of all species local populations within the regional community and synchronous dynamics among local populations of different species within the ith local community (i.e., species synchrony of the ith local community, Loreau and de Mazancourt, 2008), respectively. We can rewrite φ P→C,L (i) with d P (i, k) and s P (i, k): We define the first term of the right-hand side of the Equation S4 as the dominant-species local species synchrony of the ith local community : Then, using above description, we define the dominant-species local species synchrony of the regional community (φ d P→C,L ), that is, an aggregation of multiple local communities: Referring to the definition of local community CV, CV C,L = φ P→C,L × CV P,L Wang and Loreau, 2016;Wang and Loreau, 2014), we define the dominant-species local community CV (CV d C,L ): 1.3 Partitioning regional community synchrony into dominant and subdominant species groups The regional community synchrony (φ C,L→R ) is defined as the weighted-average synchronous dynamics among spatially separated local communities Wang and Loreau, 2016;Wang and Loreau, 2014): Using d P (i, k) and s P (i, k) mentioned above, we partition regional community synchrony into dominant (φ d C,L→A ), subdominant species groups (φ s C,L→A ), and synchronous dynamic between them (φ ds C,L→A ): Referring to the definition of regional community CV with the upscaling pathway of aggregating local communities (pathway I), CV C,R = φ C,L→R × CV C,L Wang and Loreau, 2016;Wang and Loreau, 2014), we define the dominant-species regional community CV with this upscaling pathway (CV d_C C,R ): 1.4 Partitioning regional population synchrony into dominant and subdominant species groups The regional population synchrony (φ P,L→R ) is defined as the weighted-average synchronous dynamics among spatially separated local populations of same species : where ω P,L→R (k) and φ P,L→R (k) are the contributions of population variance of the kth species to that of all species within the regional community and synchrony within the kth species among sites, respectively. We can rewrite φ P,L→R (k) with d P (i, k) and s P (i, k): We define the first term of the right-hand side of above equation as the regional population synchrony of the kth (dominant) species: Then, using the above description, we defined the dominant-species estimate of regional population synchrony (φ d P,L→R ): Referring to the definition of regional population CV, CV P,R = φ P,L→R × CV P,L , we define the dominant-species regional population CV (CV d P,R ): 1.5 Partitioning regional species synchrony into dominant and subdominant species groups The regional species synchrony (φ P→C,R ) is defined as the weighted-average synchronous dynamics among regional populations of different species : Here, v P,R (kl) is the covariance between k and l regional populations. We partition the regional species synchrony into dominant (φ d S→C,R ), subdominant species groups (φ s S→C,R ), and synchronous dynamic between them (φ ds S→C,R ) using introduced d P (i, k) and s P (i, k): Referring to the definition of regional community CV with the upscaling pathway of organizing regional populations (pathway II), CV C,R = φ P→C,R × CV P,R , we define the dominantspecies regional community CV with this upscaling pathway (CV d_P C,R ): 1.6 Comparing dominant-species regional community CVs estimated with two alternative upscaling pathways Based on a recent theoretical study , the regional community CV can be upscaled by aggregating local communities (CV C C,R ) or organizing regional populations (CV P C,R ): These descriptions (Equations S19 and S20) show that the regional community CVs estimated with two different upscaling pathways are equal to each other.
In the following part, we explain why the dominant-species regional community CV estimated with two different upscaling pathways are not equal to each other (CV d_C C,R for estimated via aggregating local communities, pathway I, and CV d_P C,R for estimated via organizing regional populations, pathway II). The dominant-species regional community CV estimated by aggregating local communities (CV d_ The dominant-species regional community CV estimated by organizing regional populations (CV d_P C,R ) has the following description: Because these two equations have either same terms or different terms ( in Equation S22), the dominant-species regional community CV estimated with two different upscaling pathways should be well correlated but not exactly the same. The denominators of these two different terms are the sum of local community variances or the sum of regional population variances, respectively. The numerators are the sum of variances (and covariances) of different dominant species within the same local communities and the sum of variances (and covariances) of the same dominant species across different local communities, respectively. These differences reflect that dominant-species regional community CVs estimated via aggregating local communities (CV d_C C,R ) and via organizing regional populations (CV d_P C,R ) focus on different dominant species within the same local communities and the same dominant species across different local communities, respectively. Owing to the potential difference, we separately reported them (Figure 3-figure supplement 1a and b). We also note that the different terms in Equation  S21 and Equation S22 can be same when considering all species. This is because, in this case, they become to , and both of them are equal to 1, resulting in the same regional community CV estimated with all species using different upscaling pathways.