mowing and grazing (Hautier et al. 2014;Fay et al. 2015;Zhao et al. 2024;Zhu et al. 2024). These disturbances profoundly impact grassland productivity, carbon cycling, and ecosystem resilience (Isbell et al. 2015;Song et al. 2019;Ma, Yan, et al. 2024). In particular, the temporal stability (hereafter referred to as stability), often calculated as the ratio of mean productivity to its temporal standard deviation, has become an essential ecological indicator reflecting ecosystem capacity to resist and recover from environmental perturbations (Tilman et al. 2006;Isbell et al. 2009;Hautier et al. 2015).

Previous research on grassland responses to global change has predominantly emphasized aboveground processes, particularly focusing on aboveground net primary productivity (ANPP) and its stability (Hautier et al. 2014;Ma et al. 2017;Wilcox et al. 2017;Shao et al. 2022;Su et al. 2022). In contrast, belowground net primary productivity (BNPP), despite often comprising a substantial proportion of total ecosystem productivity, remains understudied, largely due to the scarcity of long-term observational data (Fornara and Tilman 2008; Bardgett and van Der Putten 2014;Yang et al. 2022;Brown and Collins 2023;Ma, Zhang, et al. 2024). Emerging evidence suggests belowground processes, particularly root dynamics and soil resource allocation, may be essential for buffering ecosystems against environmental fluctuations (Berkelhammer et al. 2022;Ma et al. 2023;Li et al. 2024). Nonetheless, inconsistent findings across individual studies hinder global generalizations about whether BNPP stability is better maintained (i.e., less sensitive) than ANPP stability under global change drivers (Xu et al. 2022(Xu et al. , 2024;;Yang et al. 2022;Wang et al. 2023;Ma, Yan, et al. 2024;Ma, Zhang, et al. 2024).

The underlying ecological mechanisms driving potential differences in aboveground and belowground responses to global environmental changes remain poorly understood. Aboveground productivity and stability are often best captured by broadscale climatic forcing, especially temperature and climatic water-supply indices (e.g., mean annual precipitation and aridity), because these variables regulate growingseason energy balance and the variability of plant-available soil water, thereby shaping rapid canopy responses to interannual fluctuations (Knapp and Smith 2001;Ma et al. 2017;Gilbert et al. 2020;Zhang et al. 2022). In contrast, belowground productivity may be more strongly buffered by edaphic environments, as soil resource conditions (e.g., soil moisture and nutrient availability) integrate climatic inputs through infiltration, storage, and biogeochemical cycling and therefore often exhibit dampened temporal variability and slower dynamics (Xu et al. 2015(Xu et al. , 2024;;Eskelinen and Harrison 2015;Carroll et al. 2022;Yu et al. 2024). Additionally, relatively slow turnover rates and conservative functional traits of root systems could contribute to belowground stability by moderating responses to environmental disturbances (Bai et al. 2010;Wang et al. 2019;Hanisch et al. 2020;Yan et al. 2021). Furthermore, the ecological strategy involving carbon allocation to belowground components, which is quantified as the fraction of belowground net primary productivity ( f BNPP ), has emerged as a critical integrative mechanism potentially promoting BNPP stability (Song et al. 2019;Ma et al. 2021;Liu et al. 2024;Sun et al. 2024;Wu et al. 2024). Greater belowground allocation can enhance stability under resource-limited conditions. For example, under drought plants often increase belowground investment or rooting depth/distribution, enhancing root foraging capacity and water uptake, which may buffer interannual fluctuations in belowground production (Zhang et al. 2019). Nevertheless, explicit global-scale assessments evaluating whether belowground carbon allocation is linked to ecosystem stability and to the coordination between above-and belowground responses are currently limited.

To address these knowledge gaps, we conducted a globalscale meta-analysis encompassing 1513 observations from 113 studies across 85 grassland ecosystems worldwide. We systematically quantified the responses of ANPP and BNPP, their stability, and f BNPP to nine global change drivers: nitrogen addition (N), phosphorus addition (P), combined nitrogen and phosphorus addition (N + P), increased precipitation (+PRE), reduced precipitation (-PRE), elevated atmospheric CO 2 (eCO 2 ), warming (W), mowing (M), and grazing (G). Specifically, we aimed to address three fundamental questions: (1) Is BNPP stability generally less sensitive than ANPP stability to global change drivers? (2) Do ANPP and BNPP, along with their stability, exhibit coordinated or divergent responses to global change drivers? (3) If they show divergent responses, which environmental factors primarily govern the disparity in stability between ANPP and BNPP? We hypothesized that BNPP stability is better maintained (i.e., shows weaker responses) than ANPP stability under global change drivers due to stronger buffering by soil resources; that ANPP and BNPP exhibit coupled responses in productivity but divergent responses in stability, with the strength of abovebelowground linkages modulated by belowground carbon allocation; and that, although both components ultimately depend on soil water and nutrient availability, variation in ANPP stability would be more closely associated with broad-scale climatic water-supply indices (e.g., precipitation and aridity, used here as proxies for plant-available soil water), whereas variation in BNPP stability would be more closely associated with local soil resource conditions. Addressing these hypotheses provides essential theoretical advancements in our understanding of grassland ecosystem resilience under accelerating global environmental changes, thereby advancing our understanding of grassland ecosystem stability and functioning and long-term carbon cycling.

We conducted a comprehensive literature search across the Web of Science, ScienceDirect, and China National Knowledge Infrastructure (CNKI) databases to identify peer-reviewed publications reporting grassland manipulation experiments under global change scenarios (Table S1). Additional sources were retrieved through backward searches of relevant meta-analyses and reviews. The global change drivers considered included seven climate-related factors: nitrogen addition, phosphorus addition, combined nitrogen and phosphorus addition, increased precipitation, reduced precipitation, elevated atmospheric CO 2 , and warming, as well as two grass resource use-related factors: mowing and grazing. This search identified a total of 3554 peer-reviewed publications published between January 1980 and June 2024.

The metaanalysis focused on publications investigating the effects of global change drivers on grassland community productivity, stability, and carbon allocation. Each publication included in the analysis was evaluated based on the following criteria:

1. The experiment was conducted in natural field conditions rather than in artificially constructed grassland communities.

2. The publication reported at least one of the following variables: aboveground biomass or net primary productivity, belowground biomass or net primary productivity, or root biomass.

3. The experiment included both control and treatment groups, with at least 1 replication.

4. The experiment provided at least three consecutive years of data for temporal stability calculations.

5. Experiments conducted at different sites or at one site but with varying treatment intensities were treated as independent. Data from the same site and treatment, but reported in separate publications, were combined to extend observations over multiple years. Following systematic screening, 113 studies met the inclusion criteria. A complete list of these publications is provided in the Supporting Information. Raw data were extracted directly from texts, tables, or figures using WebPlotDigitizer 4.8 (https:// apps. autom eris. io/ wpd4/ ). Additionally, background information was collected for each study, including data sources (authors and publication year), study location (longitude and latitude), climate parameters (elevation, mean annual temperature (MAT), and mean annual precipitation (MAP)), grassland type (meadow, steppe, or desert), experimental duration, and intensity (e.g., fertilizer rate, percentage change in precipitation, incremental CO 2 concentration, temperature increase, mowing height, and stocking rate). Further, the aridity index was extracted from the WorldClim database (http:// www. world clim. org/ , 30 × 30 s resolution; Fick and Hijmans 2017), and soil variables, including water content (%), sand, silt, and clay proportions (%), pH, bulk density (g cm -3

), organic carbon content (g kg -1 ), and total nitrogen content (g kg -1 ), were obtained from the ISRIC-WISE dataset (https:// soilg rids. org/ , 30 × 30 s resolution; Batjes 2016). These variables were extracted for the period during which each study was conducted, based on coordinates of the study sites.

Because site-level time series of soil moisture were not consistently available across studies, MAP and aridity were treated as broad-scale climatic indices of water-supply conditions (i.e., proxies for plant-available soil water), while recognizing that both above-and belowground production ultimately respond to soil water availability. Temporal stability was calculated as the ratio of the mean to the standard deviation of aboveground or belowground biomass/productivity over the study period (Tilman et al. 2006). Carbon allocation (f BNPP ) was computed as the proportion of belowground productivity to total productivity: BNPP/(ANPP + BNPP) (Hui and Jackson 2006). From the 113 selected studies, a total of 1513 observations were extracted, encompassing 85 study sites that were primarily distributed across North America, Europe, and Asia (Figure 1). MAT across the studies ranged from -5.9°C to 16.6°C, and MAP ranged from 139 to 1400 mm.

Given that the most extracted contrasts did not report suitable standard deviations/standard errors (or sampling variances), and that observation-level sampling errors for derived metrics (e.g., temporal stability and f BNPP ) are not consistently defined or comparable across studies, we adopted an unweighted synthesis following previous ecological meta-analyses under similar data constraints (Chen et al. 2013;Li et al. 2021;Wang et al. 2024).

The natural log response ratio (lnRR) was used to estimate the effect size of global change drivers on above and belowground response variables (productivity, stability, and carbon allocation) as follows (Hedges et al. 1999):

where X e and X c represent the mean values of response variables in the experimental treatment and control groups, respectively. The lnRR was calculated separately for each control-treatment pair.

When multifactor experiments reported multiple single-factor treatments sharing a common control, each control-treatment pair was treated as a separate contrast (Lajeunesse 2011), while nonindependence among contrasts was accounted for in subsequent mixed-effects models.

For each response variable within each global change driver (and within any stratified subset where applicable), we calculated the mean effect size lnRR and expressed it as percentage change (%) for presentation:

Uncertainty was quantified for the group-mean effect size (i.e., within each response variable × driver group), not for individual observations. The standard error of the mean effect was calculated from the among-observation variance of lnRR within that group:

where V s is the sample variance of lnRR across observations in the group and n is the number of observations (Deng et al. 2016;Wang et al. 2024). The 95% confidence interval (CI) of the mean effect size was then calculated as:

If the 95% CI did not overlap zero, the response to global change drivers was considered statistically significant. We considered group differences to be supported when 95% CIs showed little/ no overlap. Because observation-level sampling variances were unavailable for most studies, we did not implement inversevariance weighting. For the same reason, and to ensure analytical transparency, we did not conduct publication-bias

diagnostics that require per-observation standard errors (e.g., funnel plots or Egger-type regressions).

To compare the magnitude and direction of ecosystem responses among variables, we used one-way ANOVA to test for differences in response magnitude across these variables. Post hoc multiple comparisons were performed using Tukey's HSD via the agricolae package in R (de Mendiburu and de Mendiburu 2019).

To examine relationships among key ecosystem response variables across global change drivers, we fitted linear mixedeffects models (LMMs) using the lme4 package in R (Bates et al. 2015), with significance testing based on the lmerTest package (Kuznetsova et al. 2017). In these relationship models, global change driver was included as a random intercept to account for variation across experimental contexts, and random slopes for focal predictors were added when model convergence allowed. For visualizing driver-specific patterns, we also conducted linear regressions within each global change driver and extracted slope coefficients, R 2 values, and p-values.

To test whether belowground carbon allocation modulates the coupling between above-and belowground responses, we incorporated f BNPP as an interaction term in extended models. Notably, this interaction analysis was implemented using a dual-model framework, in which we ran parallel models based on both response ratios and raw measured values from control and treatment groups, allowing us to assess the robustness of f BNPP 's moderating effects. This dual-model design was applied only to f BNPP interaction models and not to the full set of LMMs. To directly compare above-and belowground temporal stability, we additionally fitted LMMs to stability values from control plots and from treatment plots separately, with compartment (ANPP stability vs. BNPP stability) as a fixed effect. To account for nonindependence arising from multiple observations within the same study (and site/experiment when available), we included study (and site/experiment nested within study) as random intercepts and assessed fixed effects using lmerTest.

To explore the role of soil properties in regulating grassland responses to global change and simplify the analysis, principal component analysis (PCA) was conducted on soil variables using the prcomp function using the factoextra package (Kassambra and Mundt 2017) in R. The first two principal components captured most of the variation in soil properties: PC1 predominantly represented soil moisture and physical properties (e.g., soil water content, sand, silt, clay, and bulk density), while PC2 reflected soil fertility and chemical properties (e.g., pH, soil organic carbon, and total nitrogen; Table S2). Together, these components provided a concise representation of soil characteristics for subsequent analyses.

To investigate the drivers of grassland productivity, stability, and carbon allocation under diverse global change scenarios, we applied random forest models to assess the relative importance of environmental and experimental factors. Predictors included soil properties (represented by the two PCA-derived components: soil moisture and soil fertility), climate variables (MAT, MAP, and aridity), vegetation factors (grassland type and species richness), and experimental factors (intensity and duration). Separate random forest models were constructed for each global change driver to reveal context-specific response patterns. Models were implemented using the randomForest package in R (Breiman et al. 2018), with the number of variables randomly sampled at each split (mtry) optimized via grid search and repeated five-fold cross-validation to improve model performance. Variable importance was quantified based on the percent increase in mean squared error (%IncMSE) following random permutation, representing the relative importance (%) of each predictor in explaining response variation. In rare cases, predictors with negligible or counterproductive influence may show negative relative importance values, which reflect a reduction in model accuracy upon permutation due to random chance or multicollinearity.

To further interpret how individual predictors influenced response variables while accounting for nonlinear relationships and interactions, we generated partial dependence plots (PDPs) based on the fitted random forest models. PDPs depict the marginal effect of a focal predictor on the response variable by averaging predictions across the joint distribution of all other predictors, thereby providing an intuitive view of the modeled relationships. For continuous predictors, PDPs were constructed using the pdp package in R (Greenwell 2017), with the predictor range discretized into 20 intervals. To enhance interpretability and reduce overfitting artifacts, we applied locally smoothed trend (LOESS) smoothing (span = 0.8) to the raw PDP curves and shaded the outer 5% quantile regions to indicate extrapolated prediction zones. For categorical variables (e.g., grassland type), PDPs were visualized as bar plots showing the average predicted response for each category level. As PDPs estimate marginal effects averaged over the training data, category-level differences reflect their independent contributions while holding other variables constant. It is important to note that for categorical predictors with low representation in certain levels, the corresponding PDP estimates may exhibit reduced robustness. PDPs were generated separately for each global change driver to account for potential context-specific predictor-response relationships.

Given the diverse units used to report treatment intensity, a deviation standardization method was applied to scale the intensity of global change treatments (X E,i ) to a range of [0, 1] (Yu et al. 2024):

where X E,i * was the experimental treatment (E) intensity of observation i after unifying scaling into [0, 1], max 1≤j≤n {X E,j } and min 1≤j≤n {X E,j } indicates the maximum and minimum of experimental treatment intensity among n observations, respectively. All statistical analyses were conducted using R version 4.4.1 (R Core Team 2016).

Global change drivers exhibit divergent effects on grassland productivity, stability, and belowground carbon allocation (Figure 2). Overall, nutrient additions, increased precipitation, elevated CO 2 , and warming enhance both ANPP and BNPP. Specifically, ANPP increases by 27.8% under N addition, 7.1% under P addition, and 71.1% under N + P addition, while BNPP increases by 22.0%, 4.7%, and 61.3%, respectively (Figure 2a-c). Similarly, increased precipitation, elevated CO 2 , and warming raise ANPP by 54.2%, 23.1%, and 5.2%, and BNPP by 10.1%, 5.2%, and 7.2% (Figure 2d,f,g). In contrast, reduced precipitation leads to strong reductions in ANPP (-29.1%) but modest declines in BNPP (-3.6%) (Figure 2e). Notably, grassland management treatments induce contrasting responses of ANPP and BNPP: both mowing and grazing reduce ANPP (-8.1% and -28.2%, respectively) but increase BNPP (+16.8% and +25.1%) (Figure 2h,i).

Importantly, BNPP stability was generally less sensitive to global change drivers than ANPP stability. Specifically, ANPP stability declines significantly under N addition (-8.3%), P addition (-18.3%), N + P addition (-22.9%), reduced precipitation (-16.9%), warming (-7.9%), and grazing (-7.9%) (Figure 2). Only increased precipitation substantially improves ANPP stability (+39.7%). In contrast, BNPP stability shows weaker and less consistent responses. It declines significantly only under N addition (-5.1%), N + P addition (-7.8%), and mowing (-16.4%), while both increased and reduced precipitation enhance BNPP stability (+16.6% and +8.2%, respectively) (Figure 2). These results indicate weaker stability responses of BNPP than ANPP across drivers.

To further assess response magnitudes, we compared the absolute response ratios across all variables (Figure S1a). Among them, the standard deviation of ANPP (ANPP_SD) exhibits the greatest variation, followed by BNPP_SD, ANPP stability, and ANPP itself. In contrast, BNPP, BNPP stability, and f BNPP display substantially lower response magnitudes with narrower distributions. Notably, f BNPP exhibited the lowest variability, reflecting a conservative and highly constrained allocation strategy under global change.

In addition, we directly compared the absolute stability values of ANPP and BNPP under both ambient control and global-change treatment conditions (Figure S1b,c). Under ambient conditions, ANPP stability was numerically higher than BNPP stability (5.599 vs. 5.135), but the difference was not significant (p > 0.05). In contrast, under global-change treatments, BNPP stability showed

higher mean values than ANPP stability (5.144 vs. 4.689; p < 0.05).

Together, these patterns indicate that belowground stability tends to be maintained more strongly under experimental perturbations, consistent with the weaker responsiveness of BNPP stability to global change drivers observed across treatments.

Changes in ANPP are positively associated with those in BNPP across global grassland ecosystems (R 2 = 0.147, p < 0.001; Figure 3a), indicating generally coordinated productivity responses between aboveground and belowground components. This coordination is particularly pronounced under N addition, N + P addition, increased precipitation, warming, mowing, and grazing. In contrast, reduced precipitation displays an opposite trend, exhibiting a negative association and highlighting context-dependent divergence in productivity responses between ecosystem compartments. However, the stability responses of ANPP and BNPP show weak and generally nonsignificant associations (Figure 3b). This absence of strong linkage suggests fundamentally different ecological mechanisms regulating stability above-and belowground, despite their coordinated productivity responses.

To test whether belowground carbon allocation modulates these abovebelow relationships, we included f BNPP as an interaction term in additional analyses (Figure S2). Higher BNPP was associated with a stronger positive relationship between ANPP and BNPP (Figure S2a-c) and with a stronger relationship between ANPP stability and BNPP stability (Figure S2d-f). Together, these patterns support f BNPP as an integrative allocation index that covaries with the degree of coordination between above-and belowground responses under global change.

Random forest modeling reveals distinct patterns in how environmental covariates are associated with above-versus belowground productivity, stability, and carbon allocation (Figure 4). For ANPP and its stability, indices of climatic water supply (MAT, aridity, MAP) consistently rank among the most important predictors. By contrast, BNPP and especially BNPP stability are more strongly associated with soil conditions, with soil fertility and soil moisture emerging as the top-ranked predictors. Notably, BNPP stability shows a particularly strong association with soil fertility, consistent with soil resource buffering as a plausible correlate of enhanced belowground stability. The f BNPP exhibits integrative associations across climatic, soil, and vegetation domains, with MAT, species richness, aridity, and soil fertility identified as key predictors, suggesting that carbon allocation links aboveand belowground processes across multiple environmental gradients. Together, these results indicate that broad-scale climatic variability is more strongly reflected in ANPP and its stability (Figures S3 and S4), whereas variation in BNPP and BNPP stability is better captured by edaphic conditions (Figures S5 and S6), with carbon allocation integrating across both climatic and edaphic domains (Figure S7). Because soil moisture is not consistently reported across primary studies, climatic water-supply indices (e.g., MAP/aridity) provide broader coverage and should be interpreted as surrogates of plant-available water rather than evidence that productivity responds to precipitation per se.

Partial dependence plots (PDPs) further illustrate the nonlinear relationships between environmental drivers and ecosystem responses (Figure 5). PDPs indicate higher predicted BNPP stability at higher soil fertility and an increase in predicted BNPP stability across intermediate-to-high values of the soil moisture component. Conversely, climatic variables such as MAT and aridity consistently exhibit a negative influence on ANPP and its stability. The relationship between f BNPP and aridity follows a unimodal pattern, peaking under moderate drought conditions, consistent with peak predicted f BNPP under intermediate aridity in the fitted models.

Additional analyses in Figure S8 highlight the modulatory roles of secondary factors, including grassland vegetation type, species richness, experimental intensity, and duration. Grassland type significantly affects belowground stability, with meadow grasslands generally showing greater stability than desert and steppe grasslands. Species richness predominantly influences carbon allocation and ANPP stability through complex nonlinear interactions. Although experimental intensity and duration generally exert weaker influences, they nevertheless show detectable effects across all ecosystem response variables, particularly shaping belowground dynamics (Figure S8).

Collectively, these findings reinforce that grassland ecosystem responses to global change are mediated by multiple, partly distinct environmental pathways. Variation in ANPP stability is most strongly reflected in broad-scale climatic water supply, whereas variation in BNPP stability is more tightly linked to local soil-resource availability. Furthermore, the integrative role of belowground carbon allocation (f BNPP ), which simultaneously responds to climate, soil resources, and vegetation structure, represents a coordinated ecological strategy that modulates the coupling between aboveground and belowground responses under global environmental change.

Our global meta-analysis provides the first comprehensive evidence that grassland community BNPP stability is consistently less responsive to global change drivers than ANPP stability (Figure 2, Figure S1). Although ANPP and BNPP often respond in a coupled manner to these drivers, their stabilities are decoupled (Figure 3). This decoupling reveals that, although both ANPP and BNPP ultimately depend on plant-available soil water and nutrient supply, the variables that best capture cross-study variation in their stability differ in our synthesis. ANPP stability is most strongly associated with broad-scale climatic indices of water supply (e.g., MAP and aridity), whereas BNPP stability is more strongly associated with edaphic conditions such as soil fertility and soil moisture (Figure 4). Importantly, MAP/aridity are treated here as coarse climatic proxies for water-supply variability (i.e., surrogates of plantavailable soil water), rather than as a distinct mechanism from soil moisture. This divergence challenges the traditional ecological perspective that aboveand belowground ecosystem functions respond in a coordinated manner (Wardle et al. 2004;Bardgett et al. 2005;Kardol and Wardle 2010;Thakur et al. 2021), highlighting instead a greater buffering capacity of the belowground subsystem.

Several interrelated ecological mechanisms underpin the observed stronger maintenance under perturbation of belowground productivity. First, belowground systems benefit significantly from the inherent buffering capacity of soil resources. Soil moisture and nutrient availability exhibit slower dynamics and greater temporal persistence compared to aboveground climatic conditions (Reichstein et al. 2013;de Vries and Caruso 2016;Zhang, Dong, et al. 2023;Shi et al. 2024). For instance, soil moisture availability is moderated by infiltration capacity, retention properties, and evapotranspiration feedbacks, decoupling it from short-term atmospheric fluctuations (Xu et al. 2014(Xu et al. , 2015)). Similarly, soil nutrient pools regulated through microbial mineralization, sorption-desorption processes, and organic matter turnover offer a more stable resource base for root systems (Schimel and Bennett 2004;Chapin III et al. 2002). This soil-mediated resource stability may buffer belowground productivity against rapid environmental fluctuations; our synthesis shows patterns consistent with this buffering, for example the relatively stable BNPP response under drought compared with ANPP (Figure 2e). Second, root systems possess conservative life-history traits conferring temporal inertia. Roots typically exhibit longer lifespans, lower turnover rates, and less morphological plasticity than aboveground plant tissues such as leaves (Bai et al. 2010;Hanisch et al. 2020;Yan et al. 2021). These conservative traits minimize root vulnerability to transient disturbances and result in slower, more resilient responses to environmental variability (Freschet et al. 2021). Such stable root trait characteristics ensure continuous belowground carbon input and productivity, thus enhancing ecosystem resilience under fluctuating environmental conditions.

Third, soil microbial communities substantially contribute to belowground productivity stability. Although individual microbial cells turn over rapidly, the emergent properties of microbial communities, such as functional redundancy, metabolic buffering, and persistent soil organic mattermicrobe feedbacks, allow community-level functions to remain relatively stable over time (Li et al. 2025;Zhang et al. 2025). Such functional stability can buffer nutrient cycling and enhance root resource-use efficiency under variable environmental conditions. Microbial legacy effects and the persistence of key functional groups may therefore reinforce belowground productivity even as aboveground conditions fluctuate (Wagg et al. 2014;Zhang, Zhou, et al. 2023;Zhang et al. 2025).

Fourth, belowground carbon allocation (f BNPP ) emerges as a useful integrative indicator of how carbon is partitioned between above-and belowground components and how this partitioning relates to stability. Our results demonstrate that f BNPP exhibits lower variability than ANPP or BNPP across diverse global change drivers (Figure 2; Figure S1), and critically, acts as a moderator of aboveground-belowground functional linkages (Figure S2). Grassland ecosystems with higher f BNPP (i.e., a larger share of production occurring belowground) tend to show stronger coordination between above-and belowground responses, consistent with greater maintenance of ecosystem function under global change. Across stress treatments such as drought, mowing, and grazing, f BNPP often increases because BNPP is maintained or declines less than ANPP, so that a larger share of total production occurs belowground. In biomassremoval treatments in particular, part of this pattern necessarily reflects the direct reduction of ANPP by clipping or herbivory, rather than a purely physiological shift in carbon allocation. Nonetheless, the tendency to maintain root systems under stress likely helps sustain resource uptake capacity and supports post-disturbance regrowth. Overall, f BNPP is best interpreted as a tractable, system-level indicator of carbon partitioning that helps explain variation in the coordination of above-and belowground stability responses, rather than a standalone proxy for BNPP stability.

Finally, evolutionary and functional convergence may have favored greater resistance of belowground productivity in grasslands. These ecosystems evolved under disturbance-prone climatic conditions and herbivore pressures, potentially selecting for traits that enhance belowground resource capture, regeneration, and stability, such as deep rooting, clonal integration, and root nutrient reserves (Reich 2014;Klimešová and Herben 2015;Ma et al. 2018). This evolutionary selection process likely contributes to the widespread pattern of belowground stability observed across diverse grassland biomes, independent of climate zone or disturbance regime (Bardgett and van Der Putten 2014;Ma et al. 2018). Collectively, buffered soil resources, conservative root traits, stable microbial communities, strategic carbon allocation, and evolutionary selection processes interactively underpin grassland belowground productivity stability, providing a critical ecological advantage that supports sustained ecosystem function despite aboveground fluctuations.

Beyond these internal ecological mechanisms, our analysis further indicates that the environmental covariates that best capture cross-study variation in stability differ between aboveand belowground compartments. Variation in ANPP stability is most strongly associated with broad-scale climatic indices of water supply (e.g., MAP and aridity), which we interpret as coarse proxies for plant-available soil water where site-level soilmoisture time series are unavailable. In contrast, BNPP stability is more strongly associated with local edaphic conditions, particularly soil fertility and soil moisture, consistent with the close dependence of root production on belowground resource environments. Recognizing these scale-dependent correlates helps clarify why above-and belowground components of grassland productivity can exhibit divergent stability responses to environmental variability.

The relatively weak yet positive correlation between ANPP and BNPP in our synthesis (Figure 3a) is in line with many longterm studies that also report only modest coupling between above-and belowground production (Hui and Jackson 2006;Brown and Collins 2023). Several biophysical and methodological factors can contribute to this pattern. First, hydraulic constraints on maintaining a continuous transpiration stream in herbaceous tissues may limit the joint expansion of canopy height and rooting depth, so that above-and belowground growth do not scale linearly. Second, the effective rooting zone is often bounded by the typical infiltration depth of rainfall events; producing roots deeper than the mean wetting front provides little benefit, such that BNPP can be tightly constrained by event size and soil water redistribution, whereas ANPP is less constrained by the search for light in relatively open grassland canopies. Third, BNPP is frequently underestimated because most root cores sample only the upper 20-30 cm of soil, even though substantial root biomass and production can occur below these horizons (Jackson et al. 1996;Schenk and Jackson 2002). Fourth, ANPP and BNPP cannot be repeatedly measured at exactly the same locations over time: clipping and coring points are typically shifted among years to avoid reharvesting the same plants or disturbed soil, introducing additional spatial heterogeneity into temporal correlations. Finally, aboveground biomass is usually harvested from relatively large quadrats, whereas BNPP is estimated from smalldiameter soil cores, creating a strong mismatch in sampling area and variance structure. Together, these nonmechanistic factors can substantially weaken observed ANPP-BNPP correlations, even when carbon fluxes remain integrated along the soil-plant-atmosphere continuum.

To contextualize our estimates, we compared our ANPP and BNPP responses with values reported in previous global syntheses. Our effect sizes fall well within established ranges. Wilcox et al. (2017) showed that ANPP is more than twice as sensitive as BNPP to increased precipitation, consistent with our increases of 54.2% versus 10.1%. Under drought, both Wilcox et al. (2017) and Wu et al. (2011) documented strong ANPP declines but minimal BNPP changes, matching our reductions of -29.1% and -3.6%. Nutrientaddition effects (ANPP +27.8%; BNPP +22.0%) also align with Komatsu et al. (2019). Similarly, Smith et al. (2024) found disproportionately stronger aboveground impacts under drought. These comparisons demonstrate that our estimated responses agree with global empirical benchmarks and therefore reinforce the robustness of our synthesis.

The consistently less responsive stability of belowground productivity implies that grasslands possess an intrinsic resilience reservoir within their soil-root-microbial systems. This finding reshapes how ecosystem stability should be understood and assessed. Because most empirical and theoretical work has focused on aboveground responses, evaluations based solely on ANPP are likely to overestimate ecosystem vulnerability to global change. Incorporating belowground stability into assessments therefore provides a more complete and balanced view of grassland resilience. This perspective highlights the importance of belowground processes in buffering environmental fluctuations and sustaining long-term ecosystem functioning (Wang et al. 2023;Xu et al. 2024). Developing feasible indicators of belowground stability, such as root trait variation, soil organic matter dynamics, or root-microbe associations, will be an important step toward improving empirical assessments and predictive models of ecosystem responses to global change.

Despite the broad geographic and driver coverage of our synthesis, several limitations warrant consideration. First, our synthesis is constrained by the limited availability of longterm BNPP observations and belowground trait data, and these data are geographically uneven: a substantial proportion of available BNPP studies originates from temperate and alpine grasslands in China, resulting in an overrepresentation of relatively cool, C 3 -dominated systems. This geographic imbalance reflects global research effort rather than selection bias, but it underscores the need for more belowground observations in warm and semi-arid grasslands to fully evaluate global generality. Second, because our synthesis is unweighted, all observations contribute equally to the effect-size distribution. This choice was necessary because standard deviations or standard errors were unavailable for a large fraction of contrasts, and several key response variables ( f BNPP and stability metrics) are derived from multi-year time series for which sampling errors are not straightforward to define. Imputing variances under these conditions would have yielded highly uncertain weights. As a result, our confidence intervals are likely somewhat wider, and thus more conservative, than those from an ideal fully weighted meta-analysis, and we were unable to apply standard variance-based tests of publication bias to the full dataset. Third, our analysis focuses primarily on single global change drivers and does not explicitly address interactive effects among nutrient enrichment, altered precipitation, and warming. Future factorial experiments incorporating root functional traits, microbial communities, and soil fauna will be essential to advance mechanistic understanding and improve predictive models of ecosystem change.