Plant diversity’s positive effect on soil respiration diminishes with increasing productivity in global forests

Domain: chemistry
Published: 2026-02-21
Editorial outcome: major_revision

Abstract
Plant diversity is essential for regulating ecosystem functions, yet its global-scale relationship with soil respiration, a critical component of the carbon cycle, remains unexplored. While plant productivity modulates ecosystem processes, understanding how plant species richness regulates soil respiration across net primary production gradients is key to predicting carbon–climate feedbacks. Here, we integrate two global plant species richness datasets (tree species richness and vascular plant species richness, including woody and non-woody plants) with a global soil respiration dataset estimated using a deep learning model trained on 6355 field observations. We demonstrate that plant diversity enhances soil respiration in low- to mid-productivity forests (<1300 g C m−2 yr−1), but this effect diminishes in high-productivity forests. By stratifying forests along a net primary production gradient and controlling for climatic, soil, and vegetation covariates via structural equation modeling, we reveal a context-dependent role of biodiversity. Greater plant species richness amplifies soil respiration in resource-limited systems but contributes minimally in high-productivity forests, where abiotic factors exert stronger influence. Our findings provide insights into how biodiversity influences soil carbon fluxes, revealing its dynamic role in shaping ecosystem carbon dynamics across productivity gradients.

Introduction
Forest plant diversity is essential for sustaining ecosystem functions, making it critical to elucidate how the loss of this diversity might degrade soil and impair overall ecosystem health1–5. Recognizing plant diversity patterns in global forests is particularly important, as it provides valuable insights into how biodiversity contributes to ecosystem resilience and carbon dynamics across diverse environmental conditions6. Over recent decades, interest has grown in understanding the ecological importance of plant diversity for ecosystem functions3,5,7–9. Numerous studies have shown a strong, positive relationship between plant diversity and soil organic carbon (SOC) storage and vegetation productivity across various ecosystems3,10–13. As a result, biodiversity declines, driven by climate and land-use changes, could lead to important global carbon losses14. While substantial progress has been made in understanding the role of plant diversity in preserving SOC storage, its effect on soil respiration (Rs), a critical process in global carbon cycling, remains elusive.

Rs, the release of carbon dioxide (CO2) from the soil into the atmosphere, is the second largest carbon flux and a key factor in regulating global carbon cycling15,16. However, the relationship between plant diversity and Rs remains relatively unexplored, with no consensus emerging from the few studies available17–20. Most evidence has been derived from small-scale experiments that manipulate plant species richness, often focused on grassland ecosystems, although such experiments are increasingly being conducted in forests3,5,18,21–24. While it has been suggested that biodiversity experiments can effectively reflect real-world biodiversity–ecosystem functioning relationships25, others argue that such experiments, especially those of short duration, may not fully capture the complexity of natural ecosystems, where biodiversity is not randomly distributed and biotic and abiotic factors are inherently intertwined13. Furthermore, these experiments are often limited to unique locations and involve a comparison of monocultures with mixtures of relatively few species, impairing their generalizability10,23. A key limitation in understanding the global plant diversity–Rs relationship is the lack of consistent Rs observations at biodiversity monitoring sites. Addressing these limitations requires an investigation of the plant diversity–Rs relationship across global forests, which is essential for predicting future carbon dynamics and guiding biodiversity conservation efforts.

Belowground, the role of biodiversity in supporting key ecosystem functions may decline as ecosystems mature, due to higher plant biomass production and the accumulation of larger carbon and nutrient pools26. For example, a recent study showed that the relationship between soil biodiversity and ecosystem function is stronger in geologically younger and less productive systems27. In contrast, in older and more productive ecosystems, ecosystem functioning (including soil respiration) appears to become less reliant on soil biodiversity, suggesting that the accumulation of resources may reduce biodiversity’s marginal contribution to key soil processes. Aboveground, however, it remains unclear whether plant diversity follows similar patterns in influencing ecosystem functions. In particular, the extent to which plant diversity contributes more strongly to Rs in less productive ecosystems, where resource limitation may amplify diversity effects, has not been explicitly tested. If plant diversity functions analogously to belowground diversity, its influence on Rs may diminish in highly productive ecosystems. Therefore, it is crucial to consider productivity levels when examining the biodiversity–Rs relationship, as the impact of plant diversity may differ along gradients of net primary production (NPP).

In the present study, we investigate broad-scale empirical patterns between Rs, plant species richness, and productivity across global forests by integrating multiple global datasets. Drawing on evidence that soil biodiversity–ecosystem function (BEF)1,2,28 maybe stronger in less productive systems27, we hypothesized that plant diversity is positively associated with Rs, and that its influence is more pronounced in low-productivity forests. To test this, we first develop global Rs predictions using a deep learning approach29 based on 6,355 Rs field observations (Supplementary Figs. 1 and 2). Deep learning models have been increasingly adopted in Earth system science for their ability to capture complex, nonlinear relationships among interacting variables, without relying on linear assumptions30,31. We then overlay global datasets of plant diversity, represented by tree species richness (TSR)32 and vascular plant species richness (including woody and non-woody plants; VPSR)33, onto the predicted global Rs map along with environmental variables such as climate, vegetation, and soil properties (Methods). To explore the diversity–Rs relationship, we quantify the effects of plant species richness on observed Rs. To capture broader spatial patterns, we use partial dependence plots (PDPs) from model-predicted Rs at the global scale and across NPP gradients while controlling for environmental variables. We further quantify the unique contribution of TSR and VPSR to Rs by isolating their effect (θ) while controlling for other environmental factors using a multiple linear regression approach across standardized levels of NPP. Finally, to disentangle direct and indirect pathways linking environmental variables, plant species richness, and Rs across productivity gradients, we use a structural equation modeling (SEM) framework.

Results
R-modeled predictions
To predict global Rs, we used a multi-layer perceptron (MLP), a deep learning approach, trained on 6355 monthly Rs observations from the MGRsD dataset34 (Supplementary Fig. 1), covering the period 1982–2018 at the global scale. To ensure the model captured fine-scale spatiotemporal variability, we assembled a set of environmental variables at a resolution of 0.05° and at monthly intervals to better capture intra-annual variability, ensuring that the predicted annual mean Rs reflects underlying seasonal dynamics29. These predictors, including climate, soil, and vegetation variables (Methods), were extracted to match each Rs observation in the training dataset. After performing a variance inflation factor (VIF) test to check for multicollinearity, feature selection yielded a comprehensive set of 12 variables used as input to train the model. Monthly global Rs predictions were then generated for each month from 1982 to 2018 and were subsequently aggregated to derive a mean annual Rs map (Supplementary Fig. 2a). The model shows strong predictive performance, achieving a coefficient of determination (R2) of 0.76 and a relative root mean square error (RRMSE) of 0.38 using 10-fold cross-validation (Supplementary Fig. 2c). The model’s robustness was validated through leave-one-year-out (LOYO) (R2=0.74) and leave-one-site-out (LOSO) (R2=0.65) cross-validations (Supplementary Fig. 2d, e), demonstrating its ability to accurately predict Rs while capturing both spatial and temporal patterns.

Global patterns of the forest species richness–Rs relationship
To quantify the marginal contributions of species richness to Rs, we trained Random Forest (RF) models based on observed Rs, while controlling for the effects of environmental variables (Supplementary Table 1). The purpose of this analysis was not to generate new predictions, but to interpret how species richness contributes most to observed variation in Rs after controlling for confounding factors. Contributions of plant species richness were quantified using Shapley Additive Explanations (SHAP) values (Methods). SHAP values provide an interpretable measure of how each variable influences Rs across its observed range. For species richness, higher SHAP values indicate a positive contribution to Rs, meaning that locations with greater VPSR or TSR tend to have higher Rs. The analysis revealed a positive relationship between species richness and Rs, with SHAP values increasing with VPSR (R2=0.36, P<0.001) and TSR (R2=0.30, P<0.001) (Fig. 1, main panels), indicating that species richness consistently enhances Rs.

Further, to examine how species richness influences global patterns of predicted Rs, we included PDPs based on Rs predicted by the MLP model (Fig. 1, insets). PDPs were generated to visualize the marginal effect of VPSR and TSR on predicted Rs. We observed a generally positive influence of species richness on Rs, with a pronounced effect at lower richness levels that tapered off at higher levels. Importantly, the environmental covariates used in the PDPs were not included for variable selection or independent prediction of Rs, but for controlling for key environmental drivers, allowing the marginal effect of species richness to be better isolated at the global scale.

Plant diversity is an important driver of forest productivity7,10,35, though its influence is often weaker than that of abiotic factors and functional composition in naturally assembled communities36,37. Given that Rs tends to be higher in more productive forests38, analyzing this relationship along an NPP gradient can provide an effective framework for understanding these dynamics. To evaluate whether species richness influences Rs differently across distinct forest productivity levels, we stratified our analysis by dividing NPP into 10 quantile-based bins, with an equal number of data points (pixels) for each bin (Supplementary Fig. 3). Our findings reveal that species richness exhibits a more pronounced positive relationship with Rs in low-productivity forests, whereas this relationship flattens in high-productivity forests, indicating reduced additional gains from species diversity (Fig. 2). Similar results were found for 20 and 50 quantile-based NPP bins (Supplementary Fig. 4 and Supplementary Tables 2 and 3).

Spatial variability of the effect of species richness on Rs
To estimate θ (i.e., the effect of plant diversity on Rs), we adopted a multiple regression approach to control for other potential factors (see Methods for details). We found that θ estimates were mostly positive in regions with lower productivity for both VPSR and TSR (Fig. 3), suggesting that increased species richness could enhance Rs in low- and moderate-productivity forests (normalized NPP from 0 to 50). Conversely, forests with the highest productivity (normalized NPP>50) showed small or negative effects of species richness on Rs (Fig. 3a, b). The decline of this effect along the NPP gradient was significant for both VPSR (R2=0.48, P<0.001) and TSR (R2=0.49, P<0.001) (Fig. 3b, d). To propagate uncertainty in Rs estimates through our analysis, we conducted a supplementary sensitivity test using generalized least squares (GLS) regression. In this approach, each observation was weighted by the variance derived from the coefficient of variation (CV) of the Rs predictions (Methods). The resulting estimates of θ were highly consistent with the original unweighted estimates (R2=0.95 for VPSR and R2=0.82 for TSR, for weighted versus unweighted θ), indicating that incorporating Rs uncertainty did not substantially affect the inferred effect sizes and the main observed patterns (Supplementary Fig. 5). Furthermore, we compared the estimates of θ across different forest types and climate zones (Supplementary Fig. 6). The effect of species richness on Rs varied significantly across forest types and climate zones (Kruskal-Wallis test and post-hoc test, P<0.05) and was significantly different from zero in all groups (Wilcoxon test, P<0.05).

We validated the relationship between VPSR, TSR, and Rs using two methods (Supplementary Figs. 7 and 8). Significant positive relationships were found, particularly in low NPP regions, corroborating the initial patterns observed (Supplementary Fig. 7a, b). Additionally, we examined predicted Rs against observed plot-level VPSR and TSR. While there was a consistently positive trend across all productivity levels, the slope was steeper for plot-level data in low-NPP regions (Supplementary Fig. 7d). In contrast, some of the relationships between predicted Rs and plot-level VPSR were non-significant, notably for high and low NPP levels (Supplementary Fig. 7c).

The direct and indirect effects of environmental factors on Rs
In addition to the direct effect of species richness on Rs (i.e., its influence not mediated by other variables), environmental factors can affect Rs both directly and indirectly, by altering species richness, which in turn influences Rs. To evaluate the relative contributions of these pathways, we used a SEM approach. The direct effect of environmental variables on Rs was strong in low-productivity forests (β >0.75), weakened in mid-productivity forests, and increased again in high-productivity forests (Fig. 4c, e). The indirect effect of environmental factors via VPSR was most pronounced in low- to mid-productivity forests (with a coefficient of 0.04 at normalized NPP level 30) and gradually decreased along the NPP gradient approaching zero, while remaining relatively stable for TSR (Fig. 4a, d). Additionally, θ for VPSR peaked in mid-productivity forests (normalized NPP level 30) but declined as NPP increased, suggesting that VPSR has its strongest influence on Rs at intermediate productivity levels. In contrast, θ for TSR was highest in low-productivity forests and diminished in mid-productivity forests. The SEM output was validated by incorporating weights derived from Rs uncertainty, confirming the robustness of the findings (Methods and Supplementary Fig. 5e–h).

Methods
Soil respiration dataset
Rs observations used to train and predict global Rs were collected from the MGRsD dataset34, which provides measurements at monthly intervals. Numerous studies have predicted Rs using the SRDB observation data64,65; however, SRDB provides data at annual intervals, potentially overlooking seasonal variability in Rs driven by climate factors and thereby reducing prediction accuracy29. Specifically, when studies in the SRDB provided detailed daily or monthly Rs data, this finer-scale information was extracted and included in MGRsD. Therefore, the use of MGRsD enhances temporal granularity and can better capture seasonal variability and climate-driven dynamics of Rs. In MGRsD, a unique SiteID is assigned if an observation’s spatial location (i.e., unique latitude–longitude coordinates) or experimental treatment (e.g., warming, nitrogen addition, or precipitation manipulation) is different. As a result, a site generally corresponds to a specific geographic point, although different treatments at the same location are assigned distinct SiteIDs. Moreover, the same location may appear in multiple studies, leading to potential overlap in coordinates across SiteIDs from different sources. The spatial extent of each site typically represents a small, localized area defined by the measurement coordinates. The number of observations per site varied widely, ranging from a single monthly value to more than 80, as many sites were sampled repeatedly across multiple months and years. Monthly Rs values were usually aggregated from daily or sub-daily measurements. In many cases, only one day of measurement was available per month, but a substantial portion of the dataset includes multiple daily observations per month, which were averaged to generate monthly values. Although most SiteIDs contain measurements for all 12 months of the year, some records have temporal gaps, particularly during winter months in cold regions, due to logistical challenges in field sampling. Despite missing data in some months, MGRsD retains sites with limited temporal coverage, including those with only a single month of data, which helps expand spatial and ecosystem representation in the dataset. Despite the extensive spatial coverage of the MGRsD, it is worth noting that biases persist in certain regions, particularly in the Northern Hemisphere and temperate ecosystems29,66.

Our inclusion criteria were as follows: (1) Records lacking clear location, missing coordinates, or Rs information were not included, (2) measurements obtained through alkali absorption or soda lime methods were omitted, as they tend to underestimate soil CO2 fluxes67,68, (3) data from manipulative experiments such as nitrogen manipulation or fertilization were also excluded, and (4) potential duplicates and all other observations flagged as suspected problem or similar were excluded in the dataset. The final Rs dataset comprised 6355 observations spanning the period 1982 to 2018 (Supplementary Fig. 1).

Species-richness, arbuscular-mycorrhizal-tree, and vegetation datasets
To analyze species richness, we utilized datasets detailing TSR and VPSR. TSR data were obtained from the work by Liang et al.32, providing a high-resolution map (0.025°) of estimated local tree species richness (mean tree species richness per hectare) derived from a global forest inventory database from ~1.3 million sample plots. The authors used plot data from ~700 terrestrial ecoregions across the world, with each ecoregion encompassing 1800 sample plots on average, providing a robust and representative basis for their model. The dataset spans a wide range of forest types and climatic zones, making it one of the most comprehensive resources for studying global tree diversity. Moreover, its high spatial resolution and extensive coverage enable detailed analyses of species richness patterns and their links to forest ecosystem functioning across diverse biomes. While extensive, the authors acknowledged potential uncertainties in the dataset, which were mitigated to some extent. Measurement errors in forest inventories were addressed through strict taxonomic validation and data standardization, while uncertainties from the imputation of diversity maps were minimized using a three cross-validation approach. VPSR data were sourced from Sabatini et al.33, providing a global-scale assessment of vascular plant alpha diversity, specifically local species richness at 2.5 arc minute resolution. This dataset was developed using the sPlot database, which incorporates more than one million vegetation plots and covers natural ecosystems on all continents and in all biomes, ensuring that richness estimates are consistent and comparable globally. For our analysis, we used the subset built using a fine spatial grain of 1 hectare (for forests only), which is necessary for capturing the variability of plant biodiversity at local scales, particularly in diverse and complex forest ecosystems. Nonetheless, the VPSR dataset has certain limitations. As noted by the authors, plot coverage is uneven, with dense samples in Europe but sparse data from tropical regions such as the Amazon and equatorial Africa. Furthermore, several tropical forest records focused mainly on woody species, leaving gaps in overall plant composition. Although targeted data searches and stratified resampling were employed to address these issues, data gaps and associated uncertainties still exist. For AM plants, we used data from ref. 43, which provides global, high-resolution maps of aboveground vegetation biomass distribution by dominant mycorrhizal associations. In the study, the authors assembled the maps of mycorrhizal vegetation based on multiple published datasets, using a number of conversion factors to obtain per pixel values of mycorrhizal plant biomass fractions. We selected the AM plant biomass (expressed as the percentage of aboveground plant biomass of arbuscular mycorrhizal plant fraction), with a resolution of 10 arc minute. We used the MODIS NPP product (https://glass.hku.hk/archive/LAI/MODIS/0.05D/) to represent the forest productivity gradient. This dataset provides reliable productivity estimates at a resolution of 0.05°.

All variables were uniformly resampled to a resolution of 0.05° for VPSR and TSR analysis and ~0.16° for AM plant analysis, using bilinear interpolation. Forest areas were selected using a global forest dataset at a resolution of 0.05°69. Non-forest areas were masked, retaining global forests with canopy cover greater than 10%, following the Food and Agriculture Organization definition70.

Environmental variables and feature selection for Rs prediction
We initially selected 29 variables encompassing climatic, soil, and vegetation factors to train the Rs model. Climatic variables were derived from the WorldClim dataset (1901–2018) at a spatial resolution of 2.5 minutes, which included monthly minimum and maximum temperatures and total precipitation71 (https://worldclim.org/data/monthlywth.html, accessed May 2023). These high-resolution data were downscaled from CRU-TS 4.06 by the Climatic Research Unit at the University of East Anglia72. Additionally, evapotranspiration data from 1981 to 2018 were sourced from the Global Land Surface Satellite (GLASS) product at a spatial resolution of 0.05°73 (http://www.glass.umd.edu/index.html, accessed May 2023). Soil variables were obtained from SoilGrid V2.0, covering depths of 0–5 cm and 5–15 cm, and with a spatial resolution of 250 m74 (https://soilgrids.org/, accessed May 2023). These variables included soil bulk density (g cm−3), cation-exchange capacity (mmol(c) kg−1), clay content (g kg−1), coarse fragments (cm3 dm−3), pH in soil water, sand content (g kg−1), silt content (g kg−1), soil nitrogen (dg kg−1), soil organic carbon (dg kg−1), and soil organic carbon stock (t ha−1, depth 0–30 cm). Vegetation variables, sourced from the GLASS product, included the fraction of absorbed photosynthetically active radiation (FAPAR), leaf area index (LAI), gross primary productivity (GPP), and fractional vegetation cover (FVC). These were available at a spatial resolution of 0.05° for 1981–201875,76 (http://www.glass.umd.edu/index.html, accessed March 2022). For each Rs observation site, environmental variables were extracted based on the site’s longitude and latitude coordinates. When necessary, these variables were uniformly resampled to a resolution of 0.05° using the bilinear interpolation method. The final model period, incorporating all variables, spanned from 1982 to 2018.

Deep learning model description for Rs predictions
The MLP model, a fundamental deep learning structure, extends the capabilities of the feedforward neural network77. This model is renowned for its ability to capture complex, non-linear relationships, mapping capability, high parallelism, and global optimization78,79. Incorporating an excessive number of environmental factors in machine learning models can disrupt data integrity, increase multicollinearity among variables, lead to overfitting, and ultimately result in model inefficiency and deviation. To address this, we employed the VIF test to assess multicollinearity among the environmental variables. Environmental factors with high multicollinearity were eliminated, ensuring that only independent variables were included in the MLP model. This approach enhances the model’s predictive power and stability across different forest ecosystems80. The VIF test was conducted using the usdm package in R, with a threshold set at 7. To compute the mean RMSE loss and determine variable importance, we utilized the DALEX package in R81. During the feature selection process, all input variables were initially ranked based on their importance and included in the MLP modeling process. The least important indicators were sequentially removed to create new feature sets. For each iteration, a new MLP model was constructed, and this process continued until the optimal feature set was identified. The final set of variables selected for the model included total precipitation, minimum temperature, bulk density (at 5–15 cm depth), cation-exchange capacity (at 5–15 cm depth), coarse fragments (at 0–5 cm depth), clay content (at 0–5 cm depth), soil nitrogen (at 0–5 cm depth), soil organic carbon stock (at 0–30 cm depth), soil organic carbon (at 5–15 cm depth), sand content (at 0–5 cm depth), silt content (at 5–15 cm depth), and LAI (monthly data at www.glass.umd.edu).

Our MLP model comprises an input layer, two hidden layers, and an output layer. The input layer corresponds to the number of input variables. The first hidden layer contains 200 neurons, while the second hidden layer has 120 neurons. To prevent overfitting, we implemented dropout regularization, randomly dropping 10% of neurons in each hidden layer during training. The output layer consists of a single neuron dedicated to regression predictions. Two critical parameters for MLP algorithms are the number of epochs (iterations of model training) and the batch size (number of samples per training session). These parameters were optimized for model performance. The epoch count was determined based on the root mean square error (RMSE) from 10-fold cross-validation of the MLP model over an initial set of 500 epochs. The batch size was set to 64 to mitigate potential overfitting issues, given the uneven distribution and limited size of the data samples. To prevent overfitting, we employed an early stopping technique, which halts training when performance on a validation dataset begins to degrade. Although the initial training was set for 500 epochs, early stopping identified the optimal number of epochs as 117. The MLP model was trained using the Keras library (version 2.9.0) in R. Subsequently, predictions were generated for each grid cell and for each month at a spatial resolution of 0.05° for the period 1982–2018. Monthly Rs predictions were chosen to better capture intra-annual variability, ensuring that the predicted annual mean Rs reflects underlying seasonal dynamics29. We then aggregated monthly predictions to derive annual mean estimates for the period 1982 to 2018. While some environmental variables were spatially upscaled to coarser resolutions for analysis, Rs estimates were retained at their original prediction scale.

Evaluation of the MLP model performance for Rs predictions
The MLP model performance was evaluated using two metrics, the R2 and the RMSE. R2 quantifies the linear relationship between observed and predicted values, with values ranging from 0 to 1. A value closer to 1 indicates greater model efficiency. The equation for the calculation of R2 is presented as follows (Eq. 1):

R2=1−∑ni=1(yi−y^i)2∑ni=1(yi−y¯)2 (1)

where yi represents the observed value, y^i represents the predicted value, y¯ represents the mean of observed values, and n represents the number of observations.

RMSE measures the discrepancies between observed and predicted values, providing an average estimate of the model’s error. Lower RMSE values indicate better model performance. Here, we used the RRMSE, which represents the RMSE normalized by the average of the actual values. The equation for RRMSE is as follows (Eq. 2):

RRMSE=√1n∑ni=1(yi−y^i)2y¯ (2)

where yi represents the observed value, y^i represents the predicted value, and n represents the number of observations.

To assess the model’s ability to capture temporal and spatial patterns, we utilized two cross-validation methods: LOYO and LOSO82,83. In these methods, data from one year (or one site) are used for testing, while data from the remaining years (or sites) are used for training, ensuring comprehensive coverage across all years and locations. Additionally, we trained the MLP model 200 times using a Monte Carlo approach. Each time, we used a different random seed to initialize the model weights. For each grid cell, we calculated the CV of the predicted Rs values across these runs. CV values were predominantly below 100% (99.7% below CV of 100), indicating robust and consistent predictive performance. Only 0.0008% of terrestrial grid cells exhibited high uncertainty (CV > 200). The spatial pattern of uncertainty is shown in Supplementary Fig. 11.

SHAP and partial dependence plot analyses
To interpret the RF models, we applied SHAP analysis, a game-theoretic framework that quantifies the contribution of each predictor to the modeled variation in Rs while controlling for interactions among variables. SHAP values were calculated for all Rs observations using the fastshap package in R, specifically using the explain function and 1,000 Monte Carlo repetitions. To illustrate the influence of species richness on Rs, we generated SHAP dependence plots showing the relationship between richness and its corresponding SHAP values.

To visualize large-scale relationships between species richness and Rs, we generated PDPs based on Rs predicted by the MLP model. PDPs were computed using the pdp84 package in R, which estimates the marginal effect of a selected variable on the modeled outcome, while averaging over the influence of all other covariates85,86. Specifically, PDPs provide insights into the average marginal effects of variables across their entire range, revealing whether the relationships are linear, monotonic, or more complex. After conducting the PDP analysis at a global level, we performed PDP at different levels of NPP. We divided the NPP data into ten equal quantile-based bins. This method allows us to observe how the diversity–Rs relationship varies across distinct NPP ranges. Each bin contains N=25,486 data points.

It is important to note that we consistently included a comprehensive set of environmental covariates (e.g., climate, soil properties, and vegetation) across all modeling steps. While the MLP model was used to predict global Rs, some of the same covariates were also incorporated in the subsequent PDP analysis. These overlapping variables were not included for variable selection or independent prediction, but to control for key environmental influences, allowing a clearer interpretation of the marginal influence of species richness after adjusting for the effects of key environmental drivers.

Effects of species richness (θ) on Rs and mapping
Forests with high species richness, typically found in more humid regions with higher temperatures (e.g., tropical forests), often exhibit higher Rs due to favorable conditions rather than species richness alone. To remove the potential confounding effect from environmental covariates, we estimated the effect of species richness on Rs (θ) using a multiple regression approach87. We selected 14 predictor variables, including climate variables, soil physical properties, soil nutrients, and vegetation characteristics (Supplementary Table 1). Forests were pooled into different levels of NPP (normalized NPP gradient), and the effect was calculated for each bin. All data points within each NPP bin were used to fit regression models. Because we included 14 parameters in the models, a minimum sample size of ten times the number of parameters (at least 140 pixels) was required for accurate estimation of θ. We excluded NPP bins with fewer than 140 pixels to ensure robust analysis across the 100 NPP bins.

The multiple regression model was developed to estimate θ while accounting for the potential confounded effects from environmental covariates. The multiple regression model was structured as follows (Eq. 3):

log(Rs i) = α0 + θ × log(S^ i) + ∑pj=1 αj Xi,j + εi (3)

where εi ∼ i.i.d. N(0, σ2), Rs i is the soil respiration at location i, S^ i is the relative species richness at location i, calculated as the percentage of species richness relative to the maximum observed species richness within each NPP bin. Xi,j are environmental variables for j=1,...,p, α0 is the intercept term, αj are regression coefficients, and εi is the error term. Both Rs and S^ terms were log-transformed to satisfy the normality assumptions of the linear model.

To assess the robustness of our results to uncertainty in Rs estimates, we conducted a supplementary sensitivity analysis in which we explicitly propagated Rs uncertainties into the regression framework. For each pixel, we used the CV raster derived from the Rs prediction uncertainties (Supplementary Fig. 11) as a measure of observation-specific variance. We then fitted our multiple regression models using GLS in the nlme package, weighting each observation by its variance. This downweighted higher-uncertainty observations and propagated Rs uncertainty into the downstream models. The GLS-based parameter estimates and p-values were compared with those from the main analysis (multiple regression and SEM) (Supplementary Fig. 5). Although R2 values were slightly lower, the qualitative conclusions were unchanged, indicating that our results are robust to plausible levels of uncertainty in Rs predictions.

To test the robustness of our estimates obtained from the multiple regression approach, we also applied alternative modeling techniques, including two types of sequential regression models87,88. These methods were chosen because they are specifically designed to accurately separate or remove the confounding effects of environmental covariates on Rs, allowing us to validate the consistency and reliability of our results. The first sequential regression model (Eqs. 4, 5) removes the influence of environmental variables on species richness by regressing log-transformed species richness against the environmental covariates. The residuals from this regression represent the variation in species richness unexplained by the environment. These residuals are then used as a predictor in a second regression of log-transformed soil respiration, alongside the environmental variables. This approach isolates the effect of species richness on soil respiration independently from environmental drivers, while still controlling for environmental influences on Rs.

log(S^ i) = α0 + ∑pj=1 αj Xi,j + εi1 (4)

log(Rs i) = γ0 + θ × residual(log(S^ i)) + ∑pj=1 γj Xi,j + εi2 (5)

The second sequential regression model reverses the order by first regressing log-transformed soil respiration on environmental variables, then regressing the residuals from this model on log-transformed species richness. This approach estimates the effect of species richness on the component of Rs variation not explained by environmental covariates (Eqs. 6, 7):

log(Rs i) = α0 + ∑pj=1 αj Xi,j + εi1 (6)

residual(log(Rs i)) = θ × log(S^ i) + εi2 (7)

The estimated values of θ from the first sequential regression model were identical to those obtained from the multiple regression approach (Supplementary Fig. 12a, c). Additionally, estimates of θ derived from the second sequential regression model were also highly consistent with those from the multiple regression model (R2=0.95, P<0.001 and TSR R2=0.79, P<0.001) (Supplementary Fig. 12b, d), supporting the robustness of θ estimation.

To map the θ, we ensured the species richness effect was evaluated across different productivity levels. We mapped the values across the global landscape by assigning the corresponding θ to each grid cell based on its NPP value. Regressions that did not reach a significant level of 0.05 were not represented on the maps. The forest types were defined following the Global Land Cover89 and classified into five categories including deciduous broadleaf forests, deciduous needleleaf forests, evergreen broadleaf forests, evergreen needleleaf forests, and mixed forests. Climate types were defined following the Köppen-Geiger climate classification90, including arid, cold, polar, temperate, and tropical climates.

We validated the relationship between VPSR, TSR, and Rs using a two-way validation approach. This method helped account for potential spatial mismatches in the datasets and verify the consistency of our results. First, we selected observed Rs values and extracted the corresponding locations for predicted TSR and VPSR data. Second, we conducted a reciprocal analysis by examining observed VPSR and TSR values and extracting the corresponding predicted Rs data (Supplementary Fig. 7). For clarity, this validation was conducted across three defined NPP levels, defined as low level (13.1 to 555.9 g C m−2 yr−1), medium level (561.2 to 709.6 g C m−2 yr−1), and high level (710.9 to 1995.0 g C m−2 yr−1).

Structural equation modeling analysis
To explore the relationships between environmental variables, species richness, and Rs, we employed an SEM approach. Unlike multiple regression, which estimates marginal effects while controlling for covariates, SEM allows us to explicitly model the causal pathways and partition direct and indirect effects among correlated predictors. This is particularly useful in ecological systems where environmental variables simultaneously influence both biodiversity and ecosystem functioning. Therefore, we modeled the direct and indirect effects of multiple environmental drivers and plant diversity on Rs across varying productivity levels. We created a composite latent variable (environment), to represent the collective influence of climate, soil, and vegetation variables (including NPP), to reduce model complexity and account for shared variance among environmental predictors. To ensure normal distribution species richness variables (VPSR and TSR) and Rs were log-transformed. The SEM model involved two primary relationships87. First, we modeled the direct effect of the environmental composite variable on species richness (β2), hypothesizing that environmental conditions directly influence biodiversity levels. Second, we assessed how the environment directly affects Rs (β1), while also considering the indirect effects of species richness on Rs (β2 × θ). This allowed us to estimate both the direct influence of environmental variables and the contribution of species richness as a mediator. We then applied the SEM model to each NPP level. The SEM analysis was also complemented by weighted linear regression models, where observations were weighted by the inverse variance derived from Rs prediction uncertainty (Supplementary Fig. 5e–h).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study presents a compelling analysis of the relationship between plant diversity and soil respiration (Rs) across global forests, integrating extensive datasets and deep learning models to assess biodiversity’s role in soil carbon fluxes. The study’s findings—that biodiversity positively influences Rs in low-productivity forests but has a diminishing effect in high-productivity ecosystems—are significant for understanding ecosystem carbon dynamics. However, while the manuscript is methodologically rigorous and addresses an important ecological question, several key issues need to be addressed before it is suitable for publication in Nature Communications. Also, I am a bit concerned if novelty of this work meets nature communications standard.\n1. The study aggregates Rs predictions over time but does not fully explore potential seasonal variations. Given that Rs is highly sensitive to climatic fluctuations, temporal dynamics should be addressed more explicitly.\n2. The predictions were generated for each grid cell at a spatial resolution of 0.05. Other spatial resolutions can also be investigated to check the robustness and reliability of the model.\n3. There are some duplications in the text. Please correct these.", "reviewer_2": "General Comments\nThe study aims to determine the relationship between plant diversity and soil respiration (Rs), using global databases to obtain data and applying various statistical techniques to test this relationship. The authors evaluate two hypotheses: the diversity-function hypothesis and the effect of net primary production (NPP) on this relationship. The study finds that plant species richness enhances soil respiration in low- to mid-productivity forests but suggests that in highly productive forests, environmental factors might regulate soil respiration. While this study has the potential to make a significant contribution to current understanding of the biotic and abiotic drivers of soil respiration patterns globally, we have reservations about the use of modeled data in the primary data analyses. Specifically, the authors used predicted values of soil respiration and plant diversity, without specifying how these products were modeled. Our concern with this approach is that the predicted values of soil respiration are not observed values, likely containing a certain amount of error that is not accounted for in the primary analyses. Similarly, the uncertainty associated with plant diversity estimates (for either data set) is not accounted for. A second key issue that we found is how the authors upscale soil respiration and plant diversity data, which they do using bi-linear interpolation. While this may be a reasonable approach for environmental variables, plant diversity in particular does not scale linearly with spatial resolution. The scaling behaviour of soil respiration is not our area of expertise, but we would strongly encourage the authors to consult the relevant literature to develop a more quantitatively robust approach to upscaling soil respiration. For these reasons we would suggest that the authors focus on their analysis of empirically observed plant diversity and soil respiration (Supplementary Fig. 6).\nMinor comments\nThe organization of the Methods section is difficult to follow, jumping around between descriptions of data and data analysis. We suggest describing the data first, e.g., soil respiration, plant diversity, and environmental variables, and then moving on to describe the different analyses.\nWe would suggest adding soil moisture as an explanatory variable in the data analyses.\nDetailed Comments\nAbstract\nClarification of terminology: The study should clarify whether it refers to species richness or taxonomic diversity. Although these are basic concepts, specifying taxonomic diversity might improve the clarity of the information.\nL21 - 24, L25-26 We suggest deleting this phrase (\"Advancements in biodiversity research...\") and the other (\"making these datasets\"), both are focused on the method of global evaluations of biodiversity and its consequences rather than on the importance of addressing these questions at global scales.\nL37-38 Add in \"forest\" here, as now there are a lot of diversity experiments in forest ecosystems\nL39 and 42 But see Jochum et al. 2020 (https://www.nature.com/articles/s41559-020-1280-9)\nL43 One could argue that one of the main limitations is/was the availability of Rs data that is comparable over space and time.\nL60 It would be helpful for the sceptical reader to include a justification for using machine learning models in this context.\nL65 This is not exactly consistent with what is being tested, as species identity effects are thought to be unique to a species, and not a group of them. We think it would be more appropriate to call this a functional group effect, which has been widely evaluated in BEF experiments (see for example, Luo et al. 2024; https://nph.onlinelibrary.wiley.com/doi/full/10.1111/nph.19889).\nL71 Delete \"and younger\", as this is not something that is addressed directly in the subsquent analyses.\nL96 Delete one \"in\"\nL100 Strictly speaking, partial dependence plots visualize the marginal effect of one feature on the predicted outcome of a model.\nL107: Since trees are also vascular plants, it is necessary to specify that both trees and wood vascular plants are considered.\nL118 Yes - but rarely the most important in naturally-assembled communities (Ratcliffe et al. 2017, http://onlinelibrary.wiley.com/doi/10.1111/ele.12849/abstract; Von der Plas 2019, 10.1111/brv.12499)\nL130: The term vascular plants is used, but since trees are also vascular plants, it would be helpful to clarify whether the text refers explicitly to tree richness or wood richness (VPSR).\nL149 Change \"a two-way approach\" to \"using two methods\"\nL166 \"interval\" should be plural.\nMain Text\nFigure 5. Not clear from this figure if the SEM includes a path from the latent variable to NPP as well.\nL209 - 211. That the effect of richness on Rs is saturating is one thing, but it is hard to understand how the concept of functional similarity is applicable here as species' level contributions to ecosystem functioning aren't quantified directly. Perhaps re-word this to \"increases, suggesting that additional species contribute relatively less to ecosystem functioning at higher levels of diversity\".\nL215 We would suggest citing a couple relevant articles (Ratcliffe et al. 2017, Guerrero-Ramirez et al. 2018, 10.1038/s41559-017-0325-1)\nL221-223. Here it would be more appropriate to mention belowground functional traits specifically (maybe even give a couple of examples).\nL228 Aboveground plant growth, yes. Belowground growth (ie roots) and belowground NPP may be higher in colder, temperate regions.\nL248 Insert \"can\" or \"may\" before \"reduces\". It would be helpful to cite a study that has shown this phenomenon to be the case.\nL250 We suggest deleting \"which can lead to competitive exclusion\", as it has been infrequently documented (if at all; 10.1111/1365-2435.12345).\nL255 Again, it would be good to cite an empirical study showing this to be the case. The cited stuy here is an opinion piece.\nL279-281 It would be good to show how this heterogeneity has been accounted for in this study.\nL311-322 There is no mention of the spatial resolution that is used for the Rs values. This is fundamental for ensuring that this matches up with environmental variables, as well as species diversity. Is there any way of weighting observations that have low error (akin to a meta-analysis)?\nL 319: Ensure clarity regarding whether this refers to optimization steps for the latent variable or another aspect of the SEM model. If applicable, explicitly state that the section is referring to the optimization method.\nL 319: Clarify whether this section refers to the method used for optimizing the variables. If so, explicitly stating this would enhance clarity.\nL 319: Be cautious when discussing this part. Are you referring to the necessary transformations applied to the latent variable \"environment\"? If so, it would be useful to clarify that explicitly.\nL323 It is unclear from the description of the model if the environmental variables are included. Also, at this point, the data used for species diversity has yet to be introduced. I would suggest moving this section to after the description of the environmental variables and the species diversity data sets.\nL334-335: The text mentions ecosystem functions but does not provide specific examples. Including some examples would enhance clarity.\nSoil respiration (Rs) components: Since Rs is the sum of heterotrophic and autotrophic respiration, was the potential effect of microbial diversity and root biomass on respiration rates considered? A higher microbial diversity and root activity could contribute to variations in Rs.\nL342 Thank you for providing the spatial resolution of the predictions. At what resolution are the input data? If different than 0.05, can we be sure how species diversity and Rs scale?\nL351-352 Does this time range correspond to the measures of Rs and species diversity?\nL366 Now we see that this is specified. We suggest adding the time period of the Rs data set in the section describing those data.\nEnvironmental variables: There is no clear mention of a variable directly indicating soil moisture, which is a key factor in soil respiration. Was soil moisture included in the SEM? Consider specifying this, as soil moisture (e.g., soil water content, soil water potential) is a crucial factor for soil respiration (see Wang et al., 2014, https://bg.copernicus.org/articles/11/259/2014/; Dacal et al., 2019, https://www.nature.com/articles/s41559-018-0770-5).\nMGRsD Database: The respiration data are recorded on a monthly basis, but the structural equation model (SEM) was constructed at an annual scale. Were the respiration data averaged over the year? If so, how was this aggregation handled? More details about the temporal resolution of the variables included in the model would be helpful.\nMinimum forest area: The study does not specify the minimum forest size considered. At what minimum area is a site classified as a forest in this study? Forest size could influence both CO₂ concentration and plant diversity.\nL432 This section should be moved up before the description of the model.\nL406 More information is needed about TSR. One would expect that these are predicted values, in which environmental variables were used as predictors. If so, how can one be sure that the results of this paper are actually representing the true effect of species diversity on Rs? Also if this data is at scale of 0.025, how was it scaled up to 0.05, given the issues with upscaling diversity estimates? These same questions apply to the VPSR data set too.\nL417 This gives the impression that raw data were used in this analysis, which does not seem to be the case.\nL 429 Bi-linear interpolation may work for up/downscaling environmental variables, but this is definitely not the case for species diversity and probably not the case for Rs.\nL433-448 See previous comments about the sequencing of the sections in this part of the manuscript.\nL483-491 We suggest adding maps of uncertainty of model predictions (ie SD of predicted Rs) in the SI.\nL492 It is unclear where or how NPP is included in the SEM. It would be a good idea to include a path from the environmental drivers to NPP, as well as to species diversity.\nReference order: The order of references in the text does not follow a clear criterion. For example, in the first paragraph of the introduction, citation [3] appears in a sequence that does not seem to follow chronological or alphabetical order. What criterion was used for the citation order?", "reviewer_3": "This study uses analysis of global datasets to determine the impact of vascular and total plant species richness on soil respiration in forest ecosystems. Soil respiration is a major component of ecosystem function and carbon balance in forests, and it is important to know the influence of biodiversity on this given diversity patterns across the globe and changing biodiversity in the face of environmental threats. The work is appealing and interesting because it provides insight into fundamental relationships between forest diversity and soil respiration, while also adding layers of understanding through additional analyses that help understand some of the mechanistic and helps to relate the findings to ecological understanding and theory. I believe this ensures the work has appeal to those interested in ecosystem carbon cycling and biodiversity-function relations alike.\nThe work finds significant positive relationship between soil respiration and both vascular plant and total plant species richness. In short, greater biodiversity enhances soil respiration, and importantly a similar pattern was found across different forest types suggesting some universality of this important and fundamental relationship. There are sound theoretical reasons to predict that biodiversity increases soil respiration, but to see this demonstrated at a global scale is very important. What is also very interesting is that the work shows that the increase in soil respiration with biodiversity was much more pronounced at lower levels of biodiversity, with less effect as biodiversity increases at high biodiversity levels. This is consistent with the idea of functional similarity / redundancy.\nWhile climatic factors (e.g., temperature, evapotranspiration, and precipitation) emerged as critical predictors of soil respiration, the strong relationships between species richness and soil respiration held, even with climatic factors included in the modelling.\nFurthermore, the work found that species richness shows a greater positive relationship with soil respiration in lower compared to higher productivity forests, suggesting the biodiversity-soil respiration relationship is less important in the most productive forests. The work also showed that the effect of species richness on soil respiration is not affected by the proportion of plants that are AM. As the authors point out, this suggests that AM plants do not significantly alter the fundamental relationship between species richness and Rs, despite their importance for overall soil carbon dynamics and soil health.\nFinally, I like that the work includes SEM to understand the direct and indirect effects of environment in soil respiration i.e. directly on soil respiration and indirectly via environmental influencing species richness.\nAs the authors themselves point out, there are important aspects not investigated such as soil C/N ratios, and other plant-microbe interactions (N fixers, ectomycorrhizas). However, we should not expect an analysis to find answers to everything, so I do not think this is a serious limitation. We need to judge the manuscript on the merits of what it presents, not on all the other things we think could be done.\nA limitation of course is that such a global scale analysis is limited in its capacity to reveal the biological mechanisms driving the patterns seen, so there is an element of speculation in the Discussion. None-the-less, I believe that the speculation is appropriate, and its right that ecological explanations are provided that helps to link this global analysis to ecological understanding and theory. In turn, I believe this enhanced the broad appeal of the work.\nThe manuscript is well written.\nMinor edits:\nLine 8: The subject is not clear in “…where environmental drivers dominate”. What do they dominate and in what way?\nLn 36: Your citations are missing this study on biodiversity-soil respiration relationships: Johnson, D et al. 2008. Plant community composition, not diversity, regulates soil respiration in grasslands. Biol. Lett. 4, 345–348. https://doi.org/10.1098/rsbl.2008.0121\nLn 185-187: It would be helpful to the reader to remind us what you mean by direct and indirect effects here (it will save us having to look at the methods)", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "Review is attached"}, "author_response": "We sincerely thank the reviewer for their thoughtful and positive constructive feedback. We greatly appreciate your recognition of the study's methodological rigor and the significance of our findings for understanding global soil carbon dynamics.\nQ1. The study aggregates Rs predictions over time but does not fully explore potential seasonal variations. Given that Rs is highly sensitive to climatic fluctuations, temporal dynamics should be addressed more explicitly.\nResponse: We thank the reviewer for this important comment. While we agree that R is strongly influenced by seasonal climatic variation, our study focuses on the spatial relationship between biodiversity and annual mean R , rather than on the temporal dynamics of R itself. This choice reflects the nature of the biodiversity data used (e.g., VPSR and TSR), which are largely static, vary little to monthly or seasonal changes, and therefore incompatible with seasonal or monthly analyses.\nNonetheless, our R predictions were generated using a deep learning model trained on monthly environmental inputs. This modeling approach captures intra-annual variability, ensuring that the predicted annual mean R reflects underlying seasonal dynamics (as we aggregated monthly prediction maps). As detailed in our recently published study in Global Change Biology1, models using monthly inputs better incorporate seasonal fluctuations in climate and vegetation, leading to more robust annual R estimates. We have now cited this work in the Methods section.\nQ2. The predictions were generated for each grid cell at a spatial resolution of 0.05. Other spatial resolutions can also be investigated to check the robustness and reliability of the model.\nResponse: We chose a spatial resolution of 0.05° (~5.5 km) to balance computational efficiency with the need to retain spatial detail from the original biodiversity datasets. Notably, the VPSR data are provided at 2.5 arcminutes (~0.042°) and the TSR data at 0.025° resolution. Using a significantly coarser resolution (e.g., 0.25° or 0.5°) would require heavy spatial aggregation of these biodiversity layers, which would reduce local heterogeneity and potentially obscure important ecological patterns2. To avoid this information loss, we opted for a resolution that retains most of the fine-scale variation in biodiversity while remaining feasible for global-scale modeling.\nAdditionally, conducting the analysis described in Figure 3 at a coarser spatial resolution may not be viable. Our sequential regression approach involves 14 predictor variables, which requires a sufficiently large number of pixels per NPP bin to ensure robust statistical inference. For instance, at least 140 pixels per bin are needed to satisfy a 10:1 ratio of observations to predictors. At a 0.5° resolution, the number of available grid cells would drop sharply, potentially undermining the statistical power and stability of the regression analysis. Therefore, using finer resolution is essential to preserve both ecological detail and analytical robustness. Our robustness assessments, including leave-one-site-out and leave-one-year-out cross-validations (Supplementary Fig. 2d,e), demonstrate the model’s reliability across both spatial and temporal dimensions.\nLastly, to answer the reviewer’s concern, we have now included an uncertainty analysis of the R model. Specifically, we assessed predictive uncertainty by running the MLP model 200 times using a Monte Carlo approach and calculating the coefficient of variation (CV) across these repetitions (Supplementary Fig. 10). This analysis showed that the CV values were predominantly below 100 (99.8% below a CV of 100), indicating strong predictive performance. Only 0.0008% of terrestrial grid cells exhibited high uncertainty (CV > 200). This analysis has been added to the Supplementary Information (Supplementary Fig. 10) and described in the Method section.\nThe additional Methods now reads:\n“Additionally, we trained the MLP model 200 times using a Monte Carlo approach. Each time using a different random seed to initialize the model weights. For each grid cell, we calculated the coefficient of variation (CV) of the predicted R values across these runs. CV values were predominantly below 100% (99.7% below CV of 100), indicating robust and consistent predictive performance. Only 0.0008% of terrestrial grid cells exhibited high uncertainty (CV > 200). The spatial pattern of uncertainty is shown in Supplementary Fig. 10.”\nSupplementary Fig. 10 Coefficient of variation (%) of soil respiration predictions from 1982 to 2018. Repetitions were performed with Monte Carlo simulations (N = 200).\nQ3. There are some duplications in the text. Please correct these.\nResponse: We have carefully reviewed the manuscript and removed all redundant or duplicated statements to improve clarity and readability.\nREFERENCES\n1. Laffitte, B. et al. Timescale Matters: Finer Temporal Resolution Influences Driver Contributions to Global Soil Respiration. Glob. Change Biol. 31, e70118 (2025).\n2. Zhou, T. et al. Reducing spatial resolution increased net primary productivity prediction of terrestrial ecosystems: A Random Forest approach. Sci. Total Environ. 897, 165134 (2023)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have no more comments for this article. The authors have resolved my concerns.", "reviewer_2": "General comments\nWe appreciate the potential contribution of this paper to advancing current understanding of the context-dependent effects of biodiversity and environmental conditions on a key ecosystem function - soil respiration. However, the somewhat circular nature of the data analyses presented in the revised version of this manuscript continues to be problematic. In the first step, the authors use a MLP model to estimate annual Rs, using a large set of environmental variables. In a following the step, the estimated values of soil respiration are then analysed with Random Forest (Figs. 1 and 2), linear regression (Fig. 3), and structural equation (Fig. 4) models to illustrate context-dependence responses of soil respiration. Particularly in the cases of the Random Forest and structural equation models, these two models largely use the same environmental variables as used in the MLP. In the description of the MLP, it remains unclear whether the authors use estimates of soil respiration that effectively average out the effects of other features. Secondly, the error associated with the estimates of soil respiration are not propagated to the other data analyses (for an example, please see Rejou-Mechain et al. 2017). The latter point is critical for ensuring the robustness of the results of the down-stream analyses (RF, linear regression, structural equation model). Another approach would be to weigh estimates by their uncertainty/precision, as in meta-analysis or with Generalised Least Squares models. Lastly, it remains unclear what the added benefit of the structural equation model is to the key message of the paper.\nResponses to authors\n1) It is not reasonable to expect readers to have read authors' previous publication(s), particularly given the central role the related analyses play in this current manuscript. Rather than referring readers to these publications, please explain the methods in sufficient details such that a reader can understand what was done without needing to read another paper.\n2). Soil moisture (or soil water content) can be modeled with SAGA or downscaled with multiple techniques (including deep-learning models, e.g., ANN, WDL, etc.\n3) SPEI is calculated over multiple time lags, which were not not specified in the analysis), so it is difficult to assess whether the presented analysis is appropriate. I would suggest testing across multiple time lags.\n4) In the Methods section, please mention the uncertainty associated with the estimates of plant and tree diversity data sets.\n5) In the Methods section, please clarify that soil respiration is not upscaled, while other variables are.\nLiterature\nRejou-Mechain M, Tanguy A, Piponiot C, Chave J, Herault B. 2017. BIOMASS : an R package for estimating above-ground biomass and its uncertainty in tropical forests (S Goslee, Ed.). Methods in Ecology and Evolution 8.\nLine-by-Line Comments\nL6 insert \"a\" before \"deep\"\nL6 a distinction is made between \"tree richness\" and \"vascular plant richness\". Since trees are a subset of vascular plants, it is recommended to clarify this distinction or unify the terminology to avoid taxonomic ambiguity throughout the manuscript.\nL12 delete /re-word \"saturated by abiotic drivers\"\nL24 change \"production\" to \"productivity\"\nL81 delete \"and efficiency\"\nL88 and 103 is \"vegetation\" and \"vegetation structure\" the same variable(s)?\nL198-199 change \"formerly referred to as functional redundancy\" to \"also referred to as funciontal redundancy\"\nL241 change \"flow\" to \"flows\"\nL249-253 this implies greater resilience to global change drivers in high diversity areas\nL268 maybe change \"associations or nitrogen-fixed plants\" to \"or soil rhizobia associations\" (for consistency in the wording)", "reviewer_3": "I am very happy with the resubmission. The authors have done everything requested for improvement.\nI will not repeat myself from the first round of reviewing on why I think the manuscript is important, but I would like to take the opportunity to briefly summarise again the work’s value:\nSoil respiration is a major component of ecosystem function and carbon balance in forests so it is important to know the influence of biodiversity on this given diversity patterns across the globe and changing biodiversity in the face of environmental threats. The work is appealing and interesting because it provides insight into fundamental relationships between forest diversity and soil respiration. It finds significant positive relationship between soil respiration and both vascular plant and total plant species richness. So, greater biodiversity enhances soil respiration. This pattern is also found across different forest types suggesting universality of this fundamental relationship. To see this demonstrated at a global scale is very important. The work shows that the increase in soil respiration with biodiversity was much more pronounced at lower levels of biodiversity, which is consistent with the idea of functional similarity / redundancy. Furthermore, it is interesting that the work found that species richness shows a greater positive relationship with soil respiration in lower compared to higher productivity forests, suggesting the biodiversity-soil respiration relationship is less important in the most productive forests. Overall, these are really important and fundamental aspects of biodiversity-function relationships for a really key C cycle process (soil respiration) that are revealed at the global scale. This is important work.\nOne minor correction:\nLine 6-7: Add “a” so it reads “…observations, and a deep-learning….”", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "Please see attachment."}, "author_response": "We thank Reviewer #1 for their positive assessment and are pleased to hear that their concerns have been fully addressed.\nWe thank Reviewer #2 for their careful second review of our manuscript. Below, we provide a detailed, point-by-point response to all comments and suggestions raised.\n1. We acknowledge that reusing a consistent set of environmental covariates across modeling steps can raise concerns about potential circularity. However, we would like to clarify that each modeling step served a distinct purpose and was not used to guide or inform parameter tuning in subsequent steps. Specifically, the MLP model was used solely to generate spatial predictions of R , while the RF and regression-based models were used subsequently to interpret and quantify the relationships between biodiversity and R under different environmental contexts.\nThe environmental covariates were, indeed, reintroduced in the RF and regression models to control for known and important environmental drivers, enabling us to better isolate and interpret the marginal effects of biodiversity. In that sense, the inclusion of overlapping variables is not actually circular here, but a deliberate strategy to adjust for confounding and increase interpretability of biodiversity effects.\nThis comment was shared with one comment from Reviewer 5 in this round of review. As reviewer 5 suggested, we have added a clarification in the Methods section to reflect this distinction and to acknowledge the reviewer’s point regarding potential/perceived circularity.\n“It is important to note that we consistently included a comprehensive set of covariates (e.g., climate, soil properties, vegetation) across all modeling steps. Specifically, the MLP was used for R predictions and the RF models for marginal effect estimation of species richness on R . These covariates were incorporated primarily to adjust for confounding factors rather than to perform variable selection. While the reuse of the same covariates across models could theoretically introduce circularity, including the main environmental variables is necessary to better isolate the effect of species richness on R .”\n2. “In the description of the MLP, it remains unclear whether the authors use estimates of soil respiration that effectively average out the effects of other features.”\nThe MLP model predicts annual R using multiple environmental inputs but does not explicitly average out any single feature’s effect. Rather, it learns nonlinear interactions among all inputs to provide accurate point predictions. Subsequent steps, such as the RF partial dependence plots, were specifically chosen to marginalize over the effects of other variables, thus complementing the prediction step with interpretable estimates of individual variable importance.\n3. For “the error associated with the estimates of soil respiration are not propagated to the other data analyses.”\nWe agree with the reviewer that propagating the uncertainty of R estimates (from Supplementary Fig. 11) is important for evaluating the robustness of downstream analyses. Following this suggestion, we conducted a supplementary “sensitivity” analysis in which R uncertainties were incorporated into both the multiple regression and SEM frameworks. For the multiple regression analyses, we used generalized least squares (GLS) with weights based on the variance derived from the of the R predictions. For SEM analyses, we applied weighted linear models using the same variance-based weights (here inverse variance as weight). The weighted analyses produced 𝜃 estimates and SEM path coefficients that were highly consistent with the original unweighted results (see new Supplementary Fig. 5). Specifically, the weighted analyses yielded 𝜃 estimates and SEM path coefficients that were highly consistent with the original unweighted results (R2 = 0.95 for VPSR and R2 = 0.82 for TSR in multiple regressions), indicating that incorporating R uncertainty did not materially alter the observed patterns or effect sizes. We therefore retained the original unweighted analyses in the main text but present the weighted results in Supplementary Fig. 5 for scientific transparency and validation.\nWe have revised the Results and Methods sections accordingly.\nPlease see in Results:\n“To propagate uncertainty in R estimates through our analysis, we conducted a supplementary sensitivity test using generalized least squares (GLS) regressions. In this approach, each observation was weighted by the variance derived from the coefficient of variation (CV) of the R predictions (Methods). The resulting estimates of 𝜃 were highly consistent with the original unweighted estimates (R2 = 0.95 for VPSR and R2 = 0.82 for TSR, for weighted versus unweighted 𝜃), indicating that incorporating R uncertainty did not substantially affect the inferred effect sizes and the main observed patterns (Supplementary Fig. 5).”\nAnd:\n“The SEM output was validated by incorporating weights derived from R uncertainty, confirming the robustness of the findings (Methods and Supplementary Fig. 5e–h).”\nIn the Methods:\n“To assess the robustness of our results to uncertainty in R estimates, we conducted a supplementary sensitivity analysis in which we explicitly propagated R uncertainties into the regression framework. For each pixel, we used the CV raster derived from the R prediction uncertainties (Supplementary Fig. 11) as a measure of observation-specific variance. We then fitted our multiple regression models using GLS in the nlme package, weighting each observation by its variance. This downweighted higher-uncertainty observations and propagated R uncertainty into the downstream models. The GLS-based parameter estimates and p-values were compared with those from the main analysis (multiple regression and SEM) (Supplementary Fig. 5). Although R2 values were slightly lower, the qualitative conclusions were unchanged, indicating that our results are robust to plausible levels of uncertainty in R predictions.”\nImportant to note, the RF model, by contrast, was used for exploratory purposes, to generate PDPs and visualize the marginal effects of TSR and VPSR. As RF is a non-parametric, ensemble-based machine learning method that does not estimate interpretable parameters in the classical statistical sense, applying observation-level uncertainty weights is not standard practice and would not be meaningful in this context.\n4. “Lastly, it remains unclear what the added benefit of the structural equation model is to the key message of the paper.”\nWe thank the reviewer for this comment. This point was previously raised by Reviewer 5 in the previous round of review, and we addressed it in our earlier response. As noted, the SEM complements the regression analyses by allowing us to separate direct and indirect effects. Specifically, it quantifies how environmental factors influence R both directly and via species richness. This cannot be achieved with the multiple or sequential regressions alone. So, SEM can enable the use of a latent environmental variable, reduce dimensionality, and improve model clarity.\nWe believe this approach strengthens our conclusions and is clearly explain in:\nThe results section:\n“In addition to the direct effect of species richness on R (i.e., its influence not mediated by other variables), environmental factors can affect R both directly and indirectly, by altering species richness, which in turn influences R . To evaluate the relative contributions of these pathways, we used a SEM approach. The direct effect of environmental variables on R was strong in low-productivity forests (𝛽 > 0.75), weakened in mid-productivity forests, and increased again in high-productivity forests (Fig. 4c, e). The indirect effect of environmental factors via VPSR was most pronounced in low- to mid-productivity forests (with a coefficient = 0.04 at normalized NPP level 30) and gradually decreased along the NPP gradient approaching zero, while remaining relatively stable for TSR (Fig. 4a, d). Additionally, 𝜃 for VPSR peaked in mid-productivity forests (normalized NPP level 30) but declined as NPP increased, suggesting that VPSR has its strongest influence on R at intermediate productivity levels. In contrast, 𝜃 for TSR was highest in low-productivity forests and diminished in mid-productivity forests.”\nThe method section:\n“To explore the relationships between environmental variables, species richness, and R , we employed an SEM approach. Unlike multiple regression, which estimates marginal effects while controlling for covariates, SEM allows us to explicitly model the causal pathways and partition direct and indirect effects among correlated predictors. This is particularly useful in ecological systems where environmental variables simultaneously influence both biodiversity and ecosystem functioning. Therefore, we modeled the direct and indirect effects of multiple environmental drivers and plant diversity on R across varying productivity levels. We created a composite latent variable, “environment”, to represent the collective influence of climate, soil, and vegetation variables (including NPP), to reduce model complexity and account for shared variance among environmental predictors. To ensure normal distribution species richness variables (VPSR and TSR) and R were log-transformed. The SEM model involved two primary relationships34. First, we modeled the direct effect of the environmental composite variable on species richness (𝛽 ), hypothesizing that environmental conditions directly influence biodiversity levels. Second, we assessed how the environment directly affects R (𝛽 ), while also considering the indirect effects of species richness on R (𝛽 ×𝜃). This allowed us to estimate both the direct influence of environmental variables and the contribution of species richness as a mediator. We then performed the SEM model across each NPP level.”\nResponses to authors\nQ2. It is not reasonable to expect readers to have read authors' previous publication(s), particularly given the central role the related analyses play in this current manuscript. Rather than referring readers to these publications, please explain the methods in sufficient details such that a reader can understand what was done without needing to read another paper.\nResponse: We agree and have carefully reviewed our Methods section to ensure clarity. We confirm that the current manuscript includes all necessary methodological details to allow readers to understand the analyses without referring to our previous publications.\nQ3. Soil moisture (or soil water content) can be modeled with SAGA or downscaled with multiple techniques (including deep-learning models, e.g., ANN, WDL, etc.\nResponse: We thank the reviewer for this suggestion. As we mentioned in our previous response, it is challenging to find high-resolution soil moisture datasets, and existing products vary significantly across regions and seasons and even show contradictory long-term trends during a same period1. There is, therefore, a strong bias using these for global-scale analysis.\nWhile SAGA offers robust tools for soil moisture modeling, applying it at the global scale for our study period (1982–2018) is impractical due to limited availability of consistent, high-resolution input data and the extensive calibration required. Additionally, incorporating SAGA-derived soil moisture could introduce added uncertainties and reduce comparability across sites and time. Instead, the MLP model leverages widely available climate and vegetation variables that effectively capture moisture-related effects on soil respiration without explicitly modeling soil moisture. For example, a recent study by Cao et al. (2025)2 examined the relationship between tree species richness and ecosystem photosynthesis and successfully captured key ecological patterns without including soil moisture as a direct predictor. In their study, they selected temperature, precipitation, as climatic factors for determining plant growth. Additionally, Yao et al. 20213 have reconstructed soil heterotrophic respiration anomaly using two different water proxies (soil moisture and precipitation), finding broadly similar patterns with minor performance differences, supporting the suitability of our approach. Therefore, we maintain that using precipitation and related predictors provides a more reliable and consistent approach for representing moisture dynamics in our global analysis.\nQ4. SPEI is calculated over multiple time lags, which were not not specified in the analysis), so it is difficult to assess whether the presented analysis is appropriate. I would suggest testing across multiple time lags.\nResponse: In our initial response, we reported results using SPEI calculated at a 1-month lag. Following the reviewer’s suggestion, we extended this analysis to include 2- and 3-month accumulation periods to capture broader drought variability. None of these time lags showed clear or significant correlations with R , and their inclusion did not improve model performance. Please see Fig. R1.\nFig. R1 Correlation between soil respiration and Standardized Precipitation-Evapotranspiration Index (SPEI). We show the relationship for 1-, 2-, and 3-month lag.\nQ5. In the Methods section, please mention the uncertainty associated with the estimates of plant and tree diversity data sets.\nResponse: We agree and have added information on potential uncertainties for estimates of both TSR and VPSR datasets.\nQ6. In the Methods section, please clarify that soil respiration is not upscaled, while other variables are.\nResponse: We have added this information.\nLiterature\nRejou-Mechain M, Tanguy A, Piponiot C, Chave J, Herault B. 2017. BIOMASS : an R package for estimating above-ground biomass and its uncertainty in tropical forests (S Goslee, Ed.). Methods in Ecology and Evolution 8.\nLine-by-Line Comments\nL6 insert \"a\" before \"deep\"\nResponse: Done as suggested\nL6 a distinction is made between \"tree richness\" and \"vascular plant richness\". Since trees are a subset of vascular plants, it is recommended to clarify this distinction or unify the terminology to avoid taxonomic ambiguity throughout the manuscript.\nResponse: We have clarified the distinction more explicitly in both the abstract and the introduction.\nL12 delete /re-word \"saturated by abiotic drivers\"\nResponse: Done as suggested. The sentence now reads: “greater plant species richness amplifies R in resource-limited systems but contributes minimally in high-productivity forests, where abiotic factors exert stronger influence.”\nL24 change \"production\" to \"productivity\"\nResponse: Done as suggested\nL81 delete \"and efficiency\"\nResponse: Done as suggested\nL88 and 103 is \"vegetation\" and \"vegetation structure\" the same variable(s)?\nResponse: The reviewer is correct. We deleted the word structure to improve clarity.\nL198-199 change \"formerly referred to as functional redundancy\" to \"also referred to as funciontal redundancy\"\nResponse: Done as suggested\nL241 change \"flow\" to \"flows\"\nResponse: Done as suggested\nL249-253 this implies greater resilience to global change drivers in high diversity areas\nResponse: We have clarified this connection in the revised manuscript as follows:\n\"When multiple species share similar functional traits, the loss or addition of species may therefore have a smaller impact on R , decreasing the ecosystem's reliance on plant diversity for regulating soil processes and potentially diminishing or even reversing the positive relationship between species richness and Rs. Thus, this implies greater resilience to global change drivers in high-diversity forest ecosystems55,56.\"\nTo support this interpretation, we have added two citations (Isbell et al. 2015; Hisano et al. 2018) and believe this addresses your point. See Ref. 4,5.\nL268 maybe change \"associations or nitrogen-fixed plants\" to \"or soil rhizobia associations\" (for consistency in the wording)\nResponse: We agree with this suggestion and have modified the wording.\nReferences\n1. Xu, X. et al. Inconsistencies in global soil moisture products and discrepancies in their relationship with vegetation productivity. Journal of Hydrology 659, 133298 (2025).\n2. Cao, R. et al. Global evidence for a positive relationship between tree species richness and ecosystem photosynthesis. Nat. Plants 1–12 (2025) doi:10.1038/s41477-025-02046-1.\n3. Yao, Y. et al. A data-driven global soil heterotrophic respiration dataset and the drivers of its inter-annual variability. Glob. Biogeochem. Cycles 35, e2020GB006918 (2021).\n4. Isbell, F. et al. Biodiversity increases the resistance of ecosystem productivity to climate extremes. Nature 526, 574–577 (2015).\n5. Hisano, M., Searle, E. B. & Chen, H. Y. H. Biodiversity as a solution to mitigate climate change impacts on the functioning of forest ecosystems. Biol. Rev. 93, 439–456 (2018).\n(Remarks on code availability)\nWe thank Reviewer #3 for the very positive and encouraging feedback. We truly appreciate the recognition and importance of our work.\nThank you for spotting this mistake. We have added the missing article “a” on Lines 6–7\nWe thank the reviewer for the thoughtful feedback throughout the review process. We have addressed the remaining comments in this revised manuscript and appreciate the opportunity to provide clarifications. We trust that the revisions now fully reflect the reviewer’s recommendations and are suitable for final consideration."}, {"round": 3, "reviewer_comments": {"reviewer_5": "The authors have done a considerable amount of work to address concerns of all of the reviewers and the manuscript is greatly improved as a result.\nI remain unconvinced about the authors argument that the random forest analysis has merely controlled for covariates other than species richness rather than being confounded by them. This is because the random forest was used to examine the effect of variables on predicted Rs (line 106 in revision) when those predicted Rs values are the result of modeling observed Rs using many of the same covariates. The authors claim in their rebuttal, relative to the use of the random forest, that \"our intent was not to perform variable selection per se or predictions of Rs. Rather, the RF model was primarily used to estimate partial dependence plots (PDPs) for TSR and VPSR, to visualize the marginal effect of species richness on Rs while holding all other covariates constant\". And yet, in the revised manuscript, we are told about the random forest results that \"several climate-related variables, including temperature, evapotranspiration, and precipitation, emerged as dominant predictors of Rs\" (lines 108-109). This seems a contradiction as it involves both variable selection and prediction. If one regresses Y on x and produces predictions of Y, then takes those predictions and asks the question of whether they are related to x, the answer should come as no surprise. This concern would have little force if the Rs used in the random forest procedure was observed Rs, but the manuscript indicates it was Rs predicted from the MLP.\nThe authors continue to refuse to acknowledge the extreme variability in observed Rs when it is related to environmental covariates. To then examine the effect of species richness on predicted Rs has thus eliminated a massive amount of variability in how actual Rs would be related to species richness. This remains a concern that has not been addressed in the revised manuscript.\nThe authors clearly do not agree with this assessment, but I still have concerns that many of the relations uncovered by the analyses of this manuscript are the result of induced relations. The authors should include a short paragraph somewhere in the manuscript that directly addresses this issue and presents their argument as to why using covariates to produce predictions and then using many of the same covariates (and a few others including species richness) to determine what is related to predicted Rs amounts only to controlling the effects of key drivers, rather than allowing those drivers to induce relations that seem much stronger than this actually the case."}, "author_response": "We thank the reviewer for the thoughtful feedback throughout the review process. We have addressed the remaining comments in this revised manuscript and appreciate the opportunity to provide clarifications. We trust that the revisions now fully reflect the reviewer’s recommendations and are suitable for final consideration.\nThank you very much for your suggestion again. As you suggested, we re-run Random Forest directly on observed R (rather than MLP-predicted R ) , and we found a positive correlation between species richness and observed R , providing the primary evidence for the positive richness effect on R (Figure 1).\nSpecifically, the main panels now show SHAP values from a Random Forest trained directly on observed R (rather than MLP-predicted R ). The Partial Dependence Plots (PDPs) using predicted R are now included as illustrative insets to visualize diversity effects at the global scale. The inclusion of environmental covariates in both steps is intended solely to control for key drivers, not to induce correlations or perform variable selection.\nWe thank the reviewer for his detailed explanation. As we described above, we re-run Random Forest based on observed R , and found a positive correlation between species richness and observed R (Figure 1). Besides, we further did a PDP analysis to detect the marginal effect of VPSR and TSR on predicted R while averaging out other variables. It should be noted that in PDP analysis, the environmental variables were controlled to avoid their confounding effects on R but not used for R predictions or variable ranking.\nSpecifically, we have made the following changes:\n1. We now present RF trained directly on observed R as the primary analysis, ensuring that species richness effects are evaluated independently of any predicted values. The RF on observed R thus produces SHAP values reflecting marginal contributions of TSR and VPSR given environmental controls, without being influenced by the MLP predictions.\n2. The analysis reveals consistent positive relationships between species richness and R (new Figure 1), despite the higher potential variability in observed R compared to predicted R . This demonstrates that the biodiversity–R pattern is robust and not an artifact of the predictive modeling framework.\n3. The PDP analysis was applied to detect the marginal effect of VPSR and TSR on predicted R . It should be noted that in PDP analysis, the environmental variables were controlled to avoid their confounding effects on R but not used for R predictions. The directionality of richness effects (from MLP-predicted R ) remains consistent to that of observed R .\n4. To avoid any possible contradiction, we have removed our previous statement that “several climate-related variables, including temperature, evapotranspiration, and precipitation, emerged as dominant predictors of R ” (lines 108–109 in the prior version). We have also deleted the bar charts indicating the relative importance of the predictors selected.\n5. Updated the method section accordingly\nThank you again for your suggestion and we agree that it should be explicitly mentioned in the manuscript. As suggested, we have added a statement in the Results and Method (revised) sections, to address this issue. Specifically, we note that the covariates were reused not for variable selection or independent prediction, but to control for well-established environmental drivers and thereby better isolate the marginal role of species richness at global scale.\nThe new Results paragraph reads as follows:\n“Global patterns of the forest species richness–R relationship\nTo quantify the marginal contributions of species richness to R we trained Random Forest (RF) models based on observed R , while controlling for the effects of environmental variables (Supplementary Table 1). The purpose of this analysis was not to generate new predictions, but to interpret how species richness contributes to observed variation in R after controlling for confounding factors. Contributions of plant species richness were quantified using Shapley Additive Explanations (SHAP) values (Methods). SHAP values provide an interpretable measure of how each variable influences R across its observed range. For species richness, higher SHAP values indicate a positive contribution to R , meaning that locations with greater VPSR or TSR tend to have higher soil respiration after controlling for environmental covariates. The analysis revealed a positive relationship between species richness and R , with SHAP values increasing with VPSR (R2 = 0.36, P < 0.01) and TSR (R2 = 0.30, P < 0.01) (Fig. 1, main panels), indicating that species richness consistently enhances R .\nFurther, to examine how species richness influences R at the global scale, we included PDPs based on R predicted by the MLP model (Fig. 1, insets). PDPs were generated to visualize the marginal effect of VPSR and TSR on predicted R . We observed a generally positive influence of species richness on R , with a pronounced effect at lower richness levels that tapered off at higher levels. Importantly, the environmental covariates used in PDPs were not included for variable selection or independent prediction of R , but for controlling for key environmental drivers, allowing the marginal effect of species richness to be better isolated at the global scale.”\nThe new Figure 1 is shown below:\nFig. 1 Marginal effects of plant species richness on soil respiration (R ). a, SHAP-derived contribution of vascular plant species richness (VPSR) to observed R (Model efficiency: R2 = 0.75). b, SHAP-derived contribution of tree species richness (TSR) to observed R (Model efficiency: R2 = 0.73). SHAP values (y-axis) indicate the marginal contribution of species richness to R , plotted against the corresponding richness values (x-axis). SHAP values were calculated for all R observations using 1000 Monte Carlo repetitions. Black lines show linear regressions with R2 and p-values. Shaded gray bands indicate 95% confidence intervals. Insets show partial dependence plots (PDPs) of R predicted by the MLP model. The red curves show the marginal effect of VPSR and TSR on R in global forests. VPSR includes all plants (i.e., herbs, shrubs, and trees), while TSR refers specifically to tree species richness. Note: 𝑹̂ does not show the predicted values, instead it shows how the values are changing with the change in VPSR or TSR, while averaging out the effect of other independent controls.\nThe revised Method reads:\nSHAP and partial dependence plot analyses\n“To interpret the RF models, we applied SHAP analysis, a game-theoretic framework that quantifies the contribution of each predictor to the modeled variation in R while controlling for interactions among variables. SHAP values were calculated for all R observations using the fastshap package in R, specifically using the “explain” function and 1,000 Monte Carlo repetitions. To illustrate the influence of species richness on R , we generated SHAP dependence plots showing the relationship between richness and its corresponding SHAP values.\nTo visualize large-scale relationships between species richness and R , we generated PDPs based on R predicted by the MLP model. PDPs were computed using the pdp84 package in R, which estimates the marginal effect of a selected variable on the modeled outcome, while averaging over the influence of all other covariates85,86. Specifically, PDPs provide insights into the average marginal effects of variables across their entire range, revealing whether the relationships are linear, monotonic, or more complex. After conducting the PDP analysis at a global level, we performed PDP at different levels of NPP. We divided the NPP data into ten equal quantile-based bins. This method allows us to observe how the diversity–R relationship varies across distinct NPP ranges. Each bin reflects a specific NPP range with an N = 25,486 data points for each bin.”"}, {"round": 4, "reviewer_comments": {"reviewer_5": "It is a considerable improvement to use observed R_s to examine the effect of species richness. I am not suggesting additional modification to the manuscript, but I note that, as regressions, the two relations depicted in Figure 1 are not at all pleasing. The demonstration that VPSR and TSR have a positive relation with R_s is established, but the form of those relations is not. A starting point might be to use a growth curve, such as a Gompertz or Freundlich model, both of which have a decreasing rate of growth as the covariate (here, species richness) grows larger. At least for VPSR this is suggested by the negative residuals when VPSR is greater than about 250. I believe this might also be in concert with the author's hypothesis that the effect of species richness is diminished for systems with high productivity."}, "author_response": "We thank the reviewer for the suggestion to fit growth curves (e.g., Gompertz or Freundlich) to capture potential diminishing effects of species richness. We note, however, that the Y-axis in Figure 1 represents SHAP values, i.e., the marginal contributions of species richness to R predicted by the RF model, rather than directly observed R . As such, these values are not raw measurements and do not correspond to a standard dependent variable suitable for parametric regression. Our goal with Figure 1 is to provide a model-agnostic, interpretable view of the marginal effect of species richness, as derived from the RF. Using a parametric curve would impose a functional form on the SHAP values, which could be misleading. The observed flattening or slight decline at very high VPSR is naturally captured by the SHAP estimates, reflecting both the data and the model’s behavior, without forcing a specific mathematical relationship."}]

Decision letter