Deeper topsoils enhance ecosystem productivity and climate resilience in arid regions, but not in humid regions

Understanding the controlling mechanisms of soil properties on ecosystem productivity is essential for sustaining productivity and increasing resilience under a changing climate. Here we investigate the control of topsoil depth (e.g., A horizons) on long‐term ecosystem productivity. We used nationwide observations (n = 2401) of topsoil depth and multiple scaled datasets of gross primary productivity (GPP) for five ecosystems (cropland, forest, grassland, pasture, shrubland) over 36 years (1986–2021) across the conterminous USA. The relationship between topsoil depth and GPP is primarily associated with water availability, which is particularly significant in arid regions under grassland, shrubland, and cropland (r = .37, .32, .15, respectively, p < .0001). For every 10 cm increase in topsoil depth, the GPP increased by 114 to 128 g C m−2 year−1 in arid regions (r = .33 and .45, p < .0001). Paired comparison of relatively shallow and deep topsoils while holding other variables (climate, vegetation, parent material, soil type) constant showed that the positive control of topsoil depth on GPP occurred primarily in cropland (0.73, confidence interval of 0.57–0.84) and shrubland (0.75, confidence interval of 0.40–0.94). The GPP difference between deep and shallow topsoils was small and not statistically significant. Despite the positive control of topsoil depth on productivity in arid regions, its contribution (coefficients: .09–.33) was similar to that of heat (coefficients: .06–.39) but less than that of water (coefficients: .07–.87). The resilience of ecosystem productivity to climate extremes varied in different ecosystems and climatic regions. Deeper topsoils increased stability and decreased the variability of GPP under climate extremes in most ecosystems, especially in shrubland and grassland. The conservation of topsoil in arid regions and improvements of soil depth representation and moisture‐retention mechanisms are critical for carbon‐sequestration ecosystem services under a changing climate. These findings and relationships should also be included in Earth system models.


| INTRODUC TI ON
Terrestrial ecosystem productivity is essential for global food security and promoting carbon sequestration (Lorenz & Lal, 2009;Schmidhuber & Tubiello, 2007), but productivity is under pressure from climate change along with the increased frequency of fire, drought, floods, frost, and decreased biodiversity (Bellard et al., 2012;Grimm et al., 2013;Isbell et al., 2015;Wu et al., 2021;Xiao et al., 2016).Understanding different factors and mechanisms that control the variability of ecosystem primary productivity provides a scientific basis to sustain its productivity and increase its resilience.Soil is the main terrestrial reservoir of carbon, nutrients, and biota and serves as the core habitat for plant growth and the thriving of global ecosystems (Blum, 2005).Productive soils not only lead to higher crop yield (Bhardwaj et al., 2011) but also to enhanced resilience to climate change (Qiao et al., 2022).Most studies focused on biological aspects or were conducted in cropland.However, soil properties and mechanisms are important for all terrestrial ecosystems.Among these properties, soil depth is important for regulating biogeochemical and hydrological cycles and ecosystem productivity (Shangguan et al., 2017).Soil formation from bedrock weathering is often accompanied by deepening soil regolith and the buildup of organic matter in the topsoil (Phillips, 2008), with formation rates controlled by climate, topography, and organisms (Jenny, 1941).
However, soils are being lost since the past century from intensified land use change, agricultural practices, and climate extremes (Brown, 1984;Montgomery, 2007), which poses threats to ecosystem productivity (Berhe et al., 2018;Larson et al., 1983;Quinton et al., 2010).Deep topsoils have been related to high soil fertility and productivity.Studies have shown that for every 10 cm of soil loss, crop productivity dropped by 4.3% due to the loss of nutrients and organic matter (Bakker et al., 2004).For every 2.5 cm of topsoil loss in the United States, crop yield dropped up to 10% for wheat and corn (Lyles, 1975), and liming, fertilization, manure application, and other best management practices cannot compensate for the yield reduction from topsoil loss (Mielke & Schepers, 1986).Most of the studies on soil erosion have focused on cropland, but the capacity of soil depth and in particular the thickness of the A horizon to support productivity at other ecosystems and its effectiveness in different climate zones has not been extensively studied.
Topsoil depth, defined here as the thickness of the A horizon, is spatially variable (Francés & Lubczynski, 2011), due to soil forming factors (e.g., climate, vegetation, topography) and human-induced erosional and management effects (Zhang et al., 2022).The vulnerability of soil to erosion and the relationship between soil erosion and production differs by soil type (Larson et al., 1983).Plant productivity is also determined by many factors such as climate and varies for different vegetation types (O' Sullivan et al., 2020), but there are limited studies on the interactive effect of climate and soil on ecosystem productivity.
Here, we investigate the control of topsoil depth on ecosystem productivity using nationwide observations of topsoil depth and datasets of gross primary productivity (GPP).We aim to address the following questions: (1) Is the relationship of topsoil depth to GPP consistent for different ecosystems (e.g., natural vs. managed ecosystems)?(2) Is this relationship affected by different climatic conditions?(3) How strong is this relationship compared to other environmental controlling factors on GPP? ( 4) To what extent is the sensitivity of ecosystem productivity to climate change affected by topsoil depth?We hypothesize that the control of topsoil depth on ecosystem productivity is more significant in natural than in managed ecosystems and in dry and cold environments.We hypothesize that deeper topsoils can increase the resilience of ecosystem productivity under climate change and climate extremes in arid regions.
The findings of this study will advance our understanding of the role of soil in ecosystem productivity.

| Topsoil dataset
The soil dataset was obtained from the National Cooperative Soil Survey (NCSS) Soil Characterization Database which contains profile descriptions and analytical data for over 30,000 pedons collected across the USA (National Cooperative Soil Survey, 2020).The sampling locations were selected to represent each mapping unit of the SSURGO map (Soil Survey Staff, 2023), and the final dataset covered half of all the mapping units; some mapping units contained more than one pedon.The samples of each pedon were collected by horizon and analyzed using standard analytical methods (Schoeneberger et al., 2012).We selected pedons that: (1) have records of longitude and latitude coordinates and are within the conterminous United States (CONUS); (2) have explicit sampling year and were sampled after 1986 at which the GPP dataset became available; (3) have data from horizons collected from the ground surface down to at least a B horizon or a C horizon if a B horizon does not exist, and do not have reporting, buried, D (the previous letter of indicating rock), vesicular, or limnic horizons, or bi-sequences or formed in human-transported materials; (4) have measurements of soil organic carbon (SOC) content and texture for A horizons; (5) have consistent land use type since 1938 belonging to cropland, forest, grassland, pasture, or shrubland (see below for details); (6) have continuous measurements of GPP from 1986 to 2021 (see below for details).In total, 2401 pedons (about 8% of the total) fitted these criteria.
We calculated the topsoil depth in each pedon by summing up the thickness of all the A horizons.Transitional horizons (e.g., AB, AE) were not included.Topsoil SOC and texture were weight-averaged for A horizons using horizon thickness as weights.Soil texture was determined by pipette method for particles smaller than 2 mm and textural classes were based on the USDA classification (Soil Survey Staff, 2014).The total carbon was determined by dry combustion method and inorganic carbon was determined using gas chromatography with HCl addition.The SOC content was estimated by subtracting the inorganic carbon from the total carbon (Soil Survey Staff, 2014).We used a generalized additive model to explore the spatial-temporal variation of topsoil depth and found that the temporal term (year) was non-significant (p > .05,data not shown here).This indicated that the topsoil depth was stable for the study period , and we did not consider a temporal change of topsoil depth in this study.

| Productivity data
Ecosystem productivity was characterized by GPP from three different datasets.Although eddy covariance flux towers provide in situ GPP data, these data are often location-specific and cannot represent the productivity of our sample locations.Instead, we used raster GPP products which were calibrated from and correlated with the flux GPP.The first GPP dataset was derived from the Landsat satellite data with 30-m spatial resolution and 16-day temporal variation since 1986 (Robinson et al., 2018).GPP was calculated using Landsat Surface Reflectance and MOD17 algorithm (Running & Zhao, 2015).The principle of the MOD17 algorithm calculated GPP from a biome-specific light use efficiency and absorbed photosynthetically active radiation, with the former optimized using eddy covariance flux tower measurements and the latter estimated from the gridded shortwave radiation (Abatzoglou, 2013) and Landsat's Normalized Difference Vegetation Index (NDVI; Robinson et al., 2018).
The dynamic biome information was provided by the National Land Cover Database (NLCD) from 1992from , 2001from , 2006from , and 2011from (Yang et al., 2018)).Smoothing and gap-filling of the Landsat's NDVI data were conducted using a climatology driven approach (Robinson et al., 2017).For details on this GPP product, readers can refer to Robinson et al. (2018).
The GPP dataset from 1986 to 2021 was extracted to pedon locations using the Nearest Neighbor method with the Google Earth Engine.As the soil profile was selected to best represent the mapping unit, we assumed that the soil thickness of the single profile is representative within the 30-m pixel of Landsat data.This means that the short-scale variation of soil is not considered in this study.
The annual accumulated GPP was calculated for each location and then averaged for 35 years to obtain the temporal mean GPP.The pedons that have missing GPP data were removed, and descriptive statistics of GPP for different ecosystems were calculated.
Given the uncertainty of the MOD17 algorithm and the Landsat data, a multi-model comparison was conducted to minimize the effect of the uncertain GPP estimates on the subsequent analysis.
Here, we included two other GPP datasets that were generated from different principles for comparison.

Solar-induced chlorophyll fluorescence (SIF) has been recently
proposed as a better proxy for GPP (Li et al., 2018).It measures the sunlight-induced photon emission from plant chlorophyll in the range from 600 to 800 nm (Baker, 2008).SIF can be retrieved from satellite observations at certain wavelengths between 600 and 800 nm.When heat dissipation occurs at high light levels, SIF is strongly correlated with photosynthesis (Baker, 2008) and has a better performance than the traditionally used vegetation index (e.g., NDVI).The GOSIF GPP dataset used here (Li & Xiao, 2019b) was derived from a global, Orbiting Carbon Observatory-2 (OCO-2)-based SIF product (GOSIF; Li & Xiao, 2019a).GOSIF consists of 0.05° and 8-day SIF estimates globally and was based on discrete SIF observations from the OCO-2, meteorological reanalysis data (photosynthetically active radiation, air temperature, and vapor pressure deficit), and MODIS-enhanced vegetation index using a machine learning method (Li & Xiao, 2019a).GOSIF GPP is based on GOSIF and GPP-SIF relationships.To estimate GPP from SIF, GPP-SIF relationships established using GPP data and OCO-2 SIF at a number of eddy covariance flux sites were used (Li & Xiao, 2019b).GOSIF GPP consists of global GPP maps with a 0.05° spatial resolution and an 8-day time step from 2000 to 2021.
FluxCom initiative provides another way to estimate GPP globally at fine spatial and temporal resolutions.It uses an ensemble of machine learning algorithms to build relationships between eddy covariance flux tower measured GPP and remote sensing satellite data (e.g., MODIS) with and without ancillary meteorological forcings (Jung et al., 2020).In this study, we used the annual GPP from the remote sensing and meteorological data-based (RS + METEO) FluxCom product with 0.5° spatial resolution.The GPP data was presented as daily GPP data (g C m −2 day −1 ) for a specific year from 1980 to 2013 and summed for annual values.
The soil profiles collected before 2001 were removed from the GOSIF GPP dataset and the soil profiles collected after 2013 were removed from the FluxCom GPP dataset, which resulted in a total of 1657 samples and 2208 samples, respectively.The temporal mean GOSIF GPP and temporal mean FluxCom GPP were calculated by averaging annual GOSIF GPP and annual FluxCom GPP for 21 years (2001-2021) and 28 years (1986-2013), respectively.

| Land use data
We used two land cover databases to ascertain that the land use of selected sample locations was consistent over a long-term period.
The United States Geological Survey (USGS) projected land use/land cover mosaics covered 1938 to 2021 with a 250-m spatial resolution (Sohl et al., 2014(Sohl et al., , 2016)).The USGS NLCD was available for 8 years: 2001, 2004, 2006, 2008, 2011, 2013, 2016, and 2019 with a 30-m spatial resolution (Yang et al., 2018).Land use types of the sample locations were extracted for every year and the sample locations

| Environmental data
The environmental data used in this study include climate variables from TerraClimate, Köppen-Geiger climate classification, topographic variables, soil orders, soil parent materials, Watershed Boundary Dataset, and irrigation types.
TerraClimate provides monthly climate data since 1958 with a 4.6-km spatial resolution (Abatzoglou et al., 2018).We chose Ter-raClimate data instead of PRISM and Daymet datasets because the Terraclimate data provides both monthly climate data (precipitation, minimum and maximum temperature, solar radiation) and climatic water balance data (actual and potential evapotranspiration).Additionally, the annual (from 1986 to 2021) and long-term averaged climate data (precipitation, minimum, and maximum temperature) were highly correlated (Pearson correlation >.95) among these three climate datasets.Monthly precipitation (pr), minimum temperature (tmmn), maximum temperature (tmmx), actual evapotranspiration (aet), potential evapotranspiration (pet), and downward surface shortwave radiation (srad) were downloaded for sample locations for every month from 1958 to 2021 using Google Earth Engine.The annual sum of pr, aet, and pet, and annual mean of tmmn, tmmx, and srad were calculated for every year and then averaged for 64 years.
The Aridity Index (AI) was calculated (Equation 1), and it represents arid and humid conditions for AI < 1 and AI > 1, respectively (Seager et al., 2018).In our dataset, there were 1461 locations in arid regions and 940 locations in humid regions.
The global map of Köppen-Geiger climate classification was downloaded from http://www.gloh2o.org/koppen/ for present day ) at a 0.0083° resolution (Beck et al., 2018).It was derived using the method described in Peel et al. (2007) with three air temperature datasets (WorldClim V1 and V2, CHELSA V1.2) and four precipitation datasets (WorldClim V1 and V2, CHELSA V1.2, and CHPclim V1) with a 0.0083° resolution, in which the CHPclim V1.2 was downscaled from 0.05° to 0.0083° resolution using bilinear interpolation (Beck et al., 2018).The Köppen-Geiger classification has a hierarchy structure with five classes at the highest level: tropical (A), arid (B), temperate (C), cold (D), and polar (E).The class B precedes other classes, and A, C, D, E classes are mutually exclusive but not with B (Beck et al., 2018).The B was identified by mean annual precipitation <10 × P threshold , in which P threshold was determined by mean annual temperature and annual precipitation pattern.The threshold to identify temperate (C) and cold (D) was 0°C for the coldest month according to Russell (1931) and 10°C for the warmest month.The Köppen-Geiger classification was extracted to the sample locations using extract function of raster package (Hijmans et al., 2013) in R. In our dataset, samples were classified as arid (n = 564), temperate (n = 748), cold (n = 1086), and polar (n = 3).
The elevation was extracted to sample locations from the USGS 3D Elevation Program 10-meter resolution Digital Elevation Model (DEM) dataset, from which the slope was calculated using Google Earth Engine (USGS).Soil order for each pedon was determined by NRCS soil scientists at the time of sampling and if missing, the gS-SURGO (30-m resolution, https://www.nrcs.usda.gov/resources/ data-and-repor ts/descr iptio n-of-gridd ed-soil-sur ve y-geogr aphic -gssur go-database) and STATSGO (1:250,000, https://catal og.data.gov/datas et/u-s-gener al-soil-map-statsgo2) maps were used to determine soil orders.Soil temperature (1:7,500,000) and moisture   (1:9,000,000) regimes maps were obtained from USDA-NRCS.The soil parent material was obtained from the Conservation Science Partners Ecologically Relevant Geomorphology map (90-m resolution) using Google Earth Engine (Soller et al., 2009;Theobald et al., 2015).The Watershed Boundary Dataset (WBD, 1:24,000scale) provided hydrologic unit (HU) data with a scale of 1:24,000 (WBD, 2022).The watershed level (HU10) was downloaded for each sample location from the Google Earth Engine.The irrigation types (irrigated or rainfed) of the cropland samples were extracted from the 2017 MODIS irrigated agricultural dataset (250-m resolution) (Brown et al., 2019).

| Statistical analysis
Four types of analysis were used to answer our research ques- We used long-term averaged GPP and climate variables for the first three analyses and annual GPP and climate variables for the last analysis.Below are the detailed explanations.
Pearson correlation coefficients (r) were calculated between topsoil depth and temporal mean GPP across the full dataset, for each ecosystem (cropland, forest, grassland, pasture, and shrubland), and in each climatic region (arid and humid regions classified by AI and arid, temperate, and cold regions classified by Köppen-Geiger classification).It was also used to explore the relationships between temporal mean GPP and other climatic and topographic variables.The relationships between temporal mean GPP and topsoil depth were fitted using simple linear regression (Equation 2) for five ecosystems and five climatic regions using lm function in R. The a represents the intercept which is the GPP when topsoil depth is zero.The b represents the slope which is the increase of GPP for every cm increase in topsoil depth.We further evaluated the topsoil depth-GPP relationship for samples with depth <75 cm (n = 2384) in five ecosystems and five climatic regions using Pearson correlation and simple linear regression, as this covered over 99% of total samples.For the final selected pairs from different watersheds, two cases may exist: (1) the deeper topsoil aligns with greater mean GPP (our hypothesis), and (2) the shallower topsoil has greater mean GPP (alternative hypothesis).Summary statistics (e.g., the proportion of pair comparisons where the site with deeper topsoil had greater mean GPP) were calculated for these two cases in each ecosystem and each climatic region.The Wald test was used to calculate the 95% confidence interval of the proportion of pair comparisons where the site with deeper topsoil had greater mean GPP (our hypothesis).If the confidence interval of the proportion covers 0.5, it indicates that the proportion of our hypothesis is not statistically significant from that of the alternative hypothesis.The Wald test was conducted using BinomCI function in DescTools package (Signorell et al., 2023) in R. We also conducted a paired t-test to compare the mean values of GPP between deep and shallow topsoils in a specific ecosystem using t.test function in R. The normality of GPP difference was tested using Shapiro-Wilk test (Shapiro & Wilk, 1965) with shapiro.
test function in R.
Since soil and environmental factors were greatly variable across watersheds, which may affect the GPP-topsoil depth relationship of the paired samples, we further developed linear mixed-effects models and multiple linear regressions using selected paired samples to account for other soil and environmental factors.In the linear mixedeffects models, we used SOC, clay content, precipitation, minimum temperature, and an interaction term of topsoil depth (both binary data-deep and shallow and numeric values were evaluated respectively) and ecosystem or climatic regions as the fixed effects to predict GPP (Equation 3).We also added the watershed as the random effect to acknowledge the variation among watersheds.The linear mixed-effects models were developed using lmer function of the lme4 package (Bates et al., 2009) in R. The significance level of the model coefficients was calculated using lmerTest (Kuznetsova et al., 2017) and afex (Singmann et al., 2015) packages in R.
To further test whether the increase in GPP is proportional to the increase in topsoil depth of the paired samples.We developed multiple linear regressions between the change in GPP (ΔGPP) and the change in topsoil depth (ΔDepth) (Equation 4).The absolute change (Equation 5) and relative change (Equation 6) were evaluated, respectively.We also added mean depth of the paired samples, SOC, clay content, precipitation, and minimum temperature as we expected these variables to have an effect.The multiple linear regressions were developed for each of the five ecosystems and five climatic regions using lm function in R.
Structural equation modeling is a multivariate regression method to examine the causal relationships between multiple variables and uses graphics to represent the complex structure (Grace, 2006).It separates the direct and indirect causes, represents partial contributions, models the latent variables and model structure (Grace, 2006), and has been used in soil ecology to solve complex casual relationships (Eisenhauer et al., 2015).We used SEMs to understand the contributions of soil and environmental variables to ecosystem productivity.In statistics, a latent variable is the one that can only be inferred from other observed variables using a mathematical model (Dodge et al., 2003).In social or natural sciences, a latent variable can be used to represent a conceptual abstract (e.g., attitude, ability) or characterize a group, and similar to observed variables, it can be used as an independent or dependent variable in models (Bollen & Hoyle, 2012).
Here, we define productivity, light, heat, water, fertility, and topsoil as six latent variables, which were inferred from other observed or measured variables, including GPP, srad, tmmn and tmmx, pr and aet, SOC and clay content, and topsoil depth, respectively.In SEM, we define productivity as a function of five latent variables (i.e., light, heat, water, fertility, and topsoil), and the interactions between the five latent variables were not investigated here (Equation 7).The SEMs were fit for five ecosystems and five climatic regions using the sem function in lavaan package (Rosseel, 2012) in R. If the full model did not converge, some observed variables were dropped until the model converged.
The model performance of SEMs was evaluated using Confirmatory Factor Index (CFI), Tucker Lewis Index (TLI), and Root Mean Square Error of Approximation (RMSEA).The CFI measures the percent decrease in model chi-square (corrected by degrees of freedom, = 2 − df) of the User Model compared to the Baseline Model and ranges from 0 to 1 (best fit) (Equation 8).The TLI measures the percent decrease in relative chi-square (χ 2 /df) of the User Model compared to the Baseline Model and when it is closer to 1 (rounded to 1 if it is greater than 1), the model is better (Equation 9).The RMSEA measures the absolute model fitting performance and when it is smaller, the model is better (Equation 10).
To evaluate the contribution of topsoil depth to climate resilience of productivity, the annual accumulative GPP, accumulative precipitation, mean minimum temperature from 1986 to 2021 (36 years) were used.Four climatic extremes from 1986 to 2021 were considered: dry, wet, hot, and cold.To identify the dry extreme for each location, the year which received the lowest precipitation was first identified from the 36 years and the GPP of this year was obtained.
Then the percent changes in GPP and precipitation to the 36-year averaged GPP and precipitation were calculated.Similarly, the wettest, the hottest, and the coldest years were identified for each location based on the highest precipitation, the highest minimum temperature, and the lowest minimum temperature, respectively, and the percent changes in GPP, precipitation, and minimum temperature were calculated.Potential lag effects from climate were not considered in our study.The relationships between percent changes in GPP at four types of climate extremes and topsoil depth were evaluated for the five climatic regions and five ecosystems.We hypothesize that if the topsoil is deeper, ecosystem productivity is more stable and hence the percent change in GPP is smaller (closer to zero).
To  S2 and S3), it confounds the effect of topsoil depth on GPP.The relationship between GPP and topsoil depth was examined for different climatic regions (Figure 2).For the same topsoil depth, humid regions had a higher GPP (mean = 1493 g C m −2 year −1 ) than arid regions (mean = 827 g C m −2 year −1 , Figure 2c; Table S4).
The correlation of topsoil depth and GPP was more evident in arid regions (AI < 1, p < 2e-16, r = .33),while in humid regions (AI > 1), the topsoil depth was slightly negatively correlated with the GPP (7) (p = .005,r = .09,Figure 2c).For every 10 cm increase in topsoil thickness, GPP increased by 114 g C m −2 year −1 in arid regions.As for Köppen-Geiger classification, the correlation of topsoil depth and GPP was not statistically significant in the temperate and cold regions (Figure 2d), but it was statistically significant in arid regions (p < 2e-16, r = .45).For every 10 cm increase in topsoil thickness, GPP increased by 128 g C m −2 year −1 in arid regions of the Köppen-Geiger classification.This may indicate that the correlation of topsoil depth and GPP was primarily associated with water availability instead of temperature, and it was stronger in dry regions.
The samples with depth <75 cm showed similar topsoil depth-GPP relationship to that of the whole dataset (Figures S4 and S5).The correlations in shrubland (r = .35vs. .32)and dry regions (r = .35vs. .33and r = .49vs. .45)were slightly stronger than that of the whole dataset.This is likely because in shallower topsoils, the changes in topsoil depth tend to be associated more with the GPP due to lack of available water and nutrients, while in deeper topsoils, such association is low.

| Association of topsoil depth with GPP at a local scale
The association of topsoil depth and GPP at a national scale is impacted by mixing effects of climate, vegetation type, topography, and soil type.To uncover the direct association, we conducted a paired comparison of GPP of relatively deep and shallow topsoils in different watersheds by holding other variables (climate, vegetation, parent material, soil type) constant (Figure 3a; Figure S6).The paired comparison contradicted our hypothesis that deeper topsoils tend to have greater GPP (Figure 3; Figure S7).The positive association of topsoil depth with GPP occurred primarily in cropland (0.73, 95% confidence interval of 0.57-0.84,29 watersheds) and shrubland (0.75, 95% confidence interval of 0.40-0.94,6 watersheds), while in forest, grassland, and pasture, over half of the shallower topsoils had higher GPP than paired deeper topsoils (Table 1).The percent increase of GPP in deeper topsoil over paired shallower topsoil was greater in shrubland (mean = 8%) followed by pasture (mean = 5%, Figure 3; Table 1).However, the paired t-test showed that the difference of GPP in deep and shallow topsoils was marginal and nonsignificant (Table 1).In the forest, the deeper topsoil was associated with a decreased average GPP change (−1.3%).At different climatic regions, the positive association of topsoil depth with GPP was marginal and not statistically significant, but it was slightly higher in arid regions (0.61, 95% confidence interval of 0.48-0.72 and 0.41-0.78)(Table S5).
The linear mixed-effects models showed that the shallow topsoil was negatively associated with GPP (coefficient = −.01), and topsoil depth as a numeric variable was positively associated with GPP (coefficient = .0018)after accounting for the other environmental factors, but such associations were not statistically significant (Table 2).In addition, the interaction of depth and forest was statistically negatively associated with GPP (coefficient = −.01)(Table 2), which indicates that deeper topsoil was associated with smaller GPP in forest.The association of GPP with ecosystem types and climatic variables (precipitation and temperature) was stronger than that with soil characteristics (Table 2).When interacted with climatic regions, the shallow topsoil was positively associated with GPP, while depth as a numeric variable was positively associated with GPP (Table S6), but such associations were not statistically significant.The precipitation was solely associated with GPP in AI classification, while both precipitation and temperature were associated with GPP in the Köppen-Geiger classification.
The absolute and relative changes in GPP did not show clear pattern with changes of topsoil depth in five ecosystems and climatic regions (Figure S8).This was also shown in the multiple linear regressions with non-statistically significant coefficients of Δ Depth (Figure S9).But such coefficients were slightly positive in most ecosystems and climatic regions except for grassland.We also observed that the absolute changes in GPP and the coefficient of Δ Depth absolute were small in shrubland, but their relative changes and the coefficient of Δ Depth relative were markedly higher (Figure 3; Figure S9).This may indicate that relative changes in topsoil depth and GPP on small absolute numbers can be more significant than they are in shrubland.The mean depth also showed statistically significant negative association with GPP change in forest.

| Association of topsoil depth and other environmental factors with GPP
To comprehensively examine the contributions of topsoil depth (topsoil) and other soil and environmental factors (light, heat, water, fertility) to ecosystem productivity, SEMs (Figure S10) were developed for each ecosystem (Figure 4) and climatic region (Figure S11).The statistically significant variables and their coefficients are shown in the SEM results (Figure 4; Figure S11).Consistent with our linear regressions, topsoil depth was positively associated with the productivity in cropland, grassland, and shrubland (dry regions, coefficients: .09-.13), but the association was not statistically significant in soils under forest and pasture (humid regions) (Figure 4).Likewise, a positive association of topsoil depth with productivity was observed in arid regions (coefficients: .17and .33),but it was negative in humid regions (coefficient = −.26)(Figure S11).Soil fertility indicated by SOC and clay content was positively associated with the productivity in cropland and shrubland (coefficients = .06),but negatively associated with the pasture and humid regions (coefficients: −.23 and −.06).In pasture, the soils with higher clay content can be more easily compacted from field traffic and form restrictive layers for root development, which may lead to lower productivity.
In humid regions, the SOC and clay content was higher toward the north (Figure S2), which was slightly opposite to the increasing GPP toward the southeast.This may indicate that other factors (e.g., climate) may play a more important role in affecting GPP than soil fertility at this scale.
Water (pr and aet) played the most important role in enhancing productivity in all the ecosystems with coefficients ranging from .07  .36).This is probably because the main limiting factor for forests (located mostly in mountainous and humid regions) is temperature.The shrubland occurred in cold and high-elevation regions (Figure S3), so increasing temperature promoted plant productivity.Light (srad) was negatively associated with the productivity in the forest, shrubland, arid, temperate, and cold regions (coefficients: −.32 to −.15), which may be associated with drought and heat stress in arid and semi-arid regions.But light was positively associated with it in humid regions (coefficient = .11),because humid regions are often radiation-limited, and an increase in light will increase annual GPP.

| Effect of topsoil depth on climate resilience of ecosystem productivity
The control of topsoil depth on the resilience of GPP was investigated under four climatic extremes in five climatic regions and five ecosystems.The distribution of ecosystems was strongly dependent on the climatic regions, in which grassland and shrubland dominated the arid regions, while forest and pasture mainly occurred in humid regions (Figure 5).In dry years, the precipitation was 20%-80% lower than the average, and it was more severe in arid regions (Figure S12).
As shown in grassland and shrubland in arid regions, the percent change in GPP was closer to zero and the productivity was more stable when the topsoil was deeper.On the contrary, when the topsoil was shallower, the GPP change was substantial (Figure 5).In some cases (e.g., pasture in humid regions), the percent change in GPP was evenly distributed with topsoil depth, which may indicate that topsoil did not affect the GPP change in climatic extremes in these regions.
Levene's test for deep and shallow topsoils was non-significant in most cases and thus did not strongly support the hypothesis that deeper topsoils had significantly smaller changes in GPP (Fig-

ure S13
).It was likely due to that threshold values were selected from mean topsoil depth (ranging from 13 to 37 cm) and may not be able to reflect the deep topsoil cases (>40 cm).The deeper topsoils had a significantly smaller variance of the percent GPP change than shallower topsoils in cropland in arid years (Figure S13) and contributed to more stable productivity in dry extremes.However, the variance was significantly larger in deeper topsoils in forest in wet years and pasture in humid regions in wet years (Figure S13).
This may suggest that forest with deeper topsoils promoted hydraulic redistribution within the deep vadose zone and fractured rocks leading to increased variation of GPP and its behavior under climate extremes (Montaldo & Oren, 2022).In some cases, the variance difference between shallow and deep topsoils was substantial (e.g., soils of pasture in arid regions, Figure S13), but Levene's test was non-significant, which was likely affected by a small sample size.
Additionally, although the percent change decreased with topsoil depth in shrubland in arid regions (Figure 5), the deeper topsoils had a larger variance than shallower topsoils (Figure S13), which was likely affected by a few extreme samples at 30-cm topsoil depths (Figure 5).
When using Köppen-Geiger classification, a similar pattern was observed in dry and wet extremes (Figures S14-S16).The percent changes in GPP decreased with deeper topsoils in cropland in cold regions, forest in temperate regions, and shrubland in arid regions, while such a pattern was less evident in grassland and pasture, and cropland in arid and temperate regions (Figure S15).As to the hot and cold extremes, the temperature increased in hot years and decreased in cold years, which was more severe in arid and cold regions (Fig-

ure S17
).Similar to that in moisture extremes, the percent changes in GPP decreased with deeper topsoils in cropland in cold regions, forest in temperate and cold regions, grassland, and shrubland in arid regions, but it was less evident in pasture (Figure S18).Similarly, Levene's test was non-significant in most cases (Figure S19).In cropland in arid regions, the coefficients were evenly distributed with topsoil depth (Figures S15 and S18), which was likely due to extensive irrigation, so that cropland water stress may be alleviated and showed no difference from rainfed fields in terms of their relationship with topsoil depth (Figure S20).

| Assessment using GOSIF and FluxCom GPP
The same analysis was conducted using the GOSIF and FluxCom GPP (Figures S21-S38).The GOSIF and FluxCom GPP showed similar spatial distribution to the Landsat GPP (Figures S21 and   S30), but the GPP of pasture was significantly higher in these two datasets than that of Landsat GPP.The relationships between GPP and topsoil depth across climatic regions remained largely the same in GOSIF and FluxCom GPP datasets, except that the GOSIF and FluxCom GPP decreased significantly with topsoil depth in cold regions (Figures S22 and S31).The percent change in GPP under climatic extremes was also evaluated using GOSIF and FluxCom GPP.For GOSIF GPP, similarly, in dry years, it decreased, while in wet years, it increased; such a change was stronger in arid regions than in humid regions (Figure S24).The percent change in GPP decreased with topsoil depth in many ecosystems, especially grassland and

Note:
The numbers in the parentheses indicate the 95% confidence intervals of the rates where deeper topsoil aligns with greater GPP in each ecosystem using the adjusted Wald method.In the Shapiro-Wilk normality test of GPP difference, yes indicates a normal distribution, while no indicates a non-normal distribution.For paired t test of GPP, n.s.indicates that the mean values of GPP between deep and shallow topsoils are not statistically significant.
shrubland (Figures S24 and S27).However, for FluxCom GPP, it decreased in dry years and increased in wet years in nearly all the locations.It was right-triangled with topsoil depth and the percent change in GPP was closer to zero with increasing topsoil depth, particularly when it was deeper than 40 cm (Figures S33 and S35).
Under temperature extremes, the percent changes in GPP can be both positive and negative and decrease with topsoil depth in grassland and shrubland (Figure S37).Abbreviation: SOC, soil organic carbon.

| Environmental controls on topsoil depth-productivity relationship
While earlier work showed positive control of topsoil depth on crop productivity, we found that these relationships varied in different ecosystems and climatic regions.By analyzing the nationwide dataset, we found that the control of topsoil depth on plant productivity was only statistically significant in drier areas, and it was not statistically significant in wetter areas.In drier regions, GPP decreased with a decreasing soil water content, but it did not change in wetter areas (Fu et al., 2022).Therefore, when plants are under water stress, deeper topsoils aligned with higher water storage can positively contribute to plant productivity.In our dataset, humid regions were distributed in the eastern CONUS and West Coast, where deeper topsoils occurred under cropland and pasture and had lower GPP than shallower soils under forest (Figure 1), and therefore a negative association of topsoil depth with productivity was observed in humid regions (Figure 2).
The control of topsoil depth on productivity was not statistically significant in temperate regions, but it was negatively associated with the GPP in cold regions when using GOSIF and FluxCom GPP (Fig- ures S22 and S31), which contradicts our hypothesis.It was assumed that deeper topsoils with generally more SOC have a higher thermal buffering capacity (Werner et al., 2020), which might be vital in cold regions to sustain plant productivity.However, in our study, cold regions were distributed mainly in the northeastern United States, Midwest, and the Rocky Mountains with predominantly forest and cropland (Figure 1), where forest had a higher GPP but shallower topsoils than cropland.Therefore, a negative relationship of topsoil depth and productivity was observed in cold regions.In temperate regions without cold stress, the relationship between topsoil depth and ecosystem productivity was not statistically significant.
The topsoil depth-GPP relationship differed for natural (e.g., forest, grassland, shrubland) and managed (e.g., cropland, pas- relation between topsoil depth and forest productivity (Figure 1).This is likely due to sufficient water and nutrient supply in the forest ecosystem and its unique root architecture and functionality.Precipitation and topsoil SOC were high under the forest (Figure S3), and hence productivity was not restricted by water or nutrient supply.Moreover, woody forest roots include primary roots going deep into the soil and fine roots spreading laterally (Danjon et al., 2013).Therefore, in these cases, topsoil depth would have little influence on forest GPP.
The paired comparison showed that the topsoil depth-GPP relationship was stronger in shrubland, while in cropland, most deeper topsoils had greater GPP, although the GPP increase was ), its association (coefficients: .09-.33) was smaller than that of water (coefficients: .07-.87) and similar to that of heat (coefficients: .06-.39).The association of topsoil depth with plant productivity was higher in arid regions.For example, in arid regions of the Köppen-Geiger classification, the coefficient of topsoil depth to productivity (.33) reached 60% of the coefficient of water (.55).Our study also found a stronger correlation between GPP and precipitation (r = .8)than that between GPP and topsoil depth (r = −.04-.37 for five ecosystems).Despite the weaker relationship, the topsoil depth played an important role in storing water and maintaining plant productivity in arid regions.

| Topsoil depth and climate resilience
Our results showed that topsoil depth not only was associated with increased long-term averaged GPP (i.e., ecosystem productivity) in some ecosystems (cropland, grassland, shrubland) and in dry regions, but it was also associated with increased resilience of ecosystem productivity to climatic variation and extremes.In arid regions, the changes in GPP under climatic extremes were more severe, while in humid regions, the GPP was more stable with smaller changes under climatic extremes.The changes in GPP were smaller in forest and pasture (wetter regions) than that in cropland, grassland, and shrubland (drier regions).These indicate that the GPP was less stable and more easily affected by climatic extremes in drier regions.However, in arid regions (especially shrubland and grassland), as topsoil was deeper, the percent change in GPP was closer to zero and the productivity was more stable.This indicates that topsoil was associated with increased climate resilience of plant productivity especially in arid regions.
Deeper topsoils aligned with higher SOC content and water storage capacity tend to increase buffering capacity to climate change and extremes and maintain ecosystem productivity.Similarly, a recent study showed that high-quality soils reduced crop yield variability and its sensitivity to climate change by over 15% (Qiao et al., 2022).Other studies have shown that biodiversity also enhanced the resilience of ecosystem productivity to climate extremes (Isbell et al., 2015) and defined it as insurance effects which included both a buffering effect and a performance-enhancing effect (Yachi & Loreau, 1999).Accordingly, deeper topsoils tend to increase both productivity and their buffering capacity to climate change.

| Implications, limitations, and prospects
Understanding the control of topsoil depth on ecosystem productivity is important for crop production, erosion control, wasteland reclamation, ecosystem restoration, and maintaining ecosystem resilience.In this study, a comprehensive analysis was conducted to study the interrelationships between topsoil depth and ecosystem and its resilience to climate extremes.Higher ecosystem productivity and low disturbance rates over the long run also contributed to greater SOC accumulation in forests and deeper topsoils in grassland.But this is a relatively slow process, especially in terms of the annual variation of GPP and its interaction with annual climate.
Soil variables are increasingly used in earth system modeling with the increasing availability, accessibility, and accuracy of national and global soil maps (Chaney et al., 2019;Poggio et al., 2021).A constant value has been used to represent soil depth in many cases in the past, but recently more global products have been available for soil depth (Pelletier et al., 2016;Shangguan et al., 2017).Topsoil depth is an important variable, as it reflects the carbon-rich and most microbially active layer of the soil and directly affects nutrient availability and ecosystem productivity.However, there is currently no available large-scale map of topsoil depth, which would be important for understanding biogeochemical cycling and earth system model-

| CON CLUS IONS
The relationship between topsoil depth and plant productivity (GPP) was investigated for different ecosystems and climatic regions using a nationwide dataset across the CONUS.A weak positive correlation was observed between the GPP and topsoil depth in soils under grassland and shrubland.The control of topsoil depth on GPP was primarily associated with water availability, which was more significant in arid regions.Forest productivity was less associated with topsoil depth due to its higher SOC content, high precipitation, and deeper root architecture and functionality.However, the pairing of deep and shallow topsoils showed a small but non-statistically significant relationship between GPP and topsoil depth.The lack of a significant relationship may be due to different sample support of soil and GPP data (i.e., point-based and raster data), or the effects of other soil and environmental factors across watersheds.Moreover, the association of GPP with topsoil depth was smaller than that with water and similar to that with heat.The topsoil depth was also related to increased stability of ecosystem productivity to climate change in arid regions and shrubland and grassland.We conclude that topsoil depth affects ecosystem productivity and its stability and resilience to climate extremes in dry regions.Such relationship does not exist in humid regions.
that had experienced land use change from 1938 to 2021 or had different land use types based on the two databases were removed from the dataset.The final dataset included 2401 pedons belonging to five land uses: cropland (n = 699), forest (n = 802), grassland (n = 324), pasture (n = 273), and shrubland (n = 303).
tions and test our hypotheses: (1) the linear regressions and Pearson correlations were calculated to explore the topsoil depth-GPP relationship in five ecosystems and five climatic regions; (2) paired comparison of relatively shallow and deep soils was conducted in five ecosystems and five climatic regions by controlling other environmental factors constant; (3) structural equation modeling (SEM) was used to evaluate the effect of topsoil depth and other essential factors (light, heat, water, fertility) and their relative contribution to GPP; (4) the effect of topsoil depth on the resilience of GPP was investigated under four climatic extremes (dry, wet, hot, and cold).
depth-GPP relationship may be confounded by many other environmental factors (e.g., climate, topography) at a national scale, to evaluate such a relationship at a local scale with other soil, topographic, and climatic factors remaining similar, we selected one pair of relatively deep and shallow topsoils for each watershed to compare their GPP values.We used watershed as the smallest spatial unit rather than other broader classifications (e.g., ecoregions), as a watershed represents the spatial movement of water (rainfall and snow melt) across the landscape and is directly related to soil erosion and deposition.As such, soils within the same watershed often have similar hydrological patterns (e.g., hydroclimatology) and the differences in plant productivity can be easily explained by the differences in soil properties when holding other environmental factors constant (see below).In each watershed, the deep and shallow topsoils were not determined by absolute depth but relative depth difference between them.The depth difference of the pair should be greater than 3 cm if they are shallow (≤15 cm) or greater than 5 cm if they are deep (>15 cm).The paired topsoils have (1) the same climate conditions (long-term averaged precipitation and temperature), land use, soil order, soil temperature and moisture regimes, parent material, and textural class; (2) similar SOC content, elevation, and slope; (3) different GPP values; (4) are within 4-km distance.If multiple pairs fulfilled the criteria for a watershed, the pair that has the most similar soil properties and the largest depth difference was selected.A total of 103 pairs were selected.

|
compare the percent change in GPP of shallow and deep topsoils, we calculated the mean topsoil depth in each climatic and ecosystem region as a threshold value and then split the data into shallow topsoils and deep topsoils.We compared the mean and variance of the percent change in GPP of shallow and deep topsoils.The mean value represents an absolute comparison of the percent GPP change, and we hypothesize that the mean percent change in GPP was closer to zero (no matter positive or negative) in deeper topsoils.The variance value represents the spread of the percent GPP change, and we hypothesize that deeper topsoils had more stable ecosystem productivity with climate change and thus the variance is smaller than that in shallower topsoils.The t-test and Levene's test were used to assess the homogeneity of mean and variance, respectively, in these two groups (shallow and deep topsoils) using t.test function and leveneTest function in car package(Fox et al., 2012) in R. If the homogeneity was not rejected in these tests, it indicates that the topsoil depth had non-significant effect on the percent change in GPP under climate extremes.Given the uncertainty and noises of the Landsat GPP data, we used another two datasets (GOSIF GPP from 2001 to 2021 and Flux-Com GPP from 1986 to 2013) to repeat two analyses: (1) the simple linear regression between temporal mean GPP and topsoil depth in five different ecosystems and climatic regions was calculated; (2) the percent changes of GPP in four climatic extremes (dry, wet, hot, and cold) were calculated, and their relationships with topsoil depth were evaluated for five climatic regions and five ecosystems.The t-test and Levene's test were also used to assess the homogeneity of mean and variance of percent GPP change for two soil depth groups (shallow and deep) in five climatic regions and five ecosystems.Topsoil depth and GPP in different ecosystems and climate zonesOn the national scale, areas with deeper topsoil were not identical to the areas with higher GPP (Figure1a,b).Deeper topsoils primarily occurred in the Midwest (Figure1a) under cropland (mean depth = 27 cm, mean GPP = 1249 g C m −2 year −1 , Figure1c,d; TableS1), while the highest GPP occurred in the east and along the West Coast (Figure1b) under forest (mean GPP = 1466 g C m −2 year −1 , Figure1c,e) with topsoils of about 15 cm (Figure1a,d).The western CONUS was dominated by shrubland and grassland (Figure1c) which had shallow topsoils (mean depth = 14 and 18 cm, Figure1a,d) and lower GPP (mean GPP = 243 and 563 g C m −2 year −1 , Figure1b,e).A weak positive correlation was observed between the GPP and topsoil depth for all data, but it varied for different ecosystems (Figure1f).A stronger positive correlation existed for grassland (r = .37,p = 7e-12) and shrubland (r = .32,p = 2e-8) with a weak correlation in cropland (r = .15,p = 4e-5), and no clear relationships for forest and pasture (r = −.04 and .02respectively, p > .05).As GPP varied significantly with climate and topography (Figures S1-S3; Tables

F
Topsoil depth and gross primary productivity (GPP) across the conterminous United States (CONUS) and in five ecosystems.(a) Topsoil depth measurements across the CONUS.(b) Landsat GPP dataset across the CONUS.(c) The distribution of sampled locations in five ecosystems.(d) The distribution of topsoil depth across five ecosystems.Numbers in parentheses indicate the sample size for each ecosystem.(e) The distribution of Landsat GPP across five ecosystems.The black dots in (d) and (e) indicate the mean values in each ecosystem.The black lines in (d) and (e) indicate the mean values across all the ecosystems.(f) The relationship between GPP and topsoil depth for the five ecosystems.The Pearson correlations between GPP and topsoil depth are .15,−.04, .37,.02,.32 for cropland, forest, grassland, pasture, and shrubland, respectively.Shadows indicate the 95% confidence intervals.
The changes in GPP due to climatic extremes were relatively small in F I G U R E 3 Paired comparison of gross primary productivity (GPP) for 103 watersheds.(a) In each watershed, one pair of relatively deep (orange dots) and shallow (blue dots) topsoils were selected to compare their GPP values.The deep and shallow topsoils were not determined by absolute depth but relative depth difference between them.The depth difference in the pair should be greater than 3 cm if they are shallower (<15 cm) or greater than 5 cm if they are deeper (>15 cm).The paired topsoils have (1) the same climate conditions, ecosystem type, soil order, soil temperature and moisture regimes, parent material, and soil textural class; (2) similar organic carbon content, elevation, and slope; (3) different GPP values; (4) are within 4-km distance.Each pair is connected by a gray line, in which the line width indicates the depth difference in the pair.If an orange dot is on top of a blue dot, it indicates the deeper topsoil aligns with greater GPP; otherwise, a shallower topsoil aligns with greater GPP.The pairs from different watersheds in each ecosystem were ranked from left to right in the x-axis by their mean annual precipitation.(b) The distribution of absolute change in GPP between deep and shallow topsoils in each watershed.(c) The distribution of relative change in GPP between deep and shallow topsoils in each watershed.The colors of the dots in (b) and (c) indicate five ecosystems, while the sizes of the dots indicate the GPP difference between deep and shallow topsoils.The two maps in (b) and (c) represent the two scenarios: (1) the deeper topsoil aligns with greater GPP (GPP differences are positive); (2) a shallower topsoil aligns with greater GPP (GPP differences are negative).
ture) ecosystems.Previous studies and topsoil removal experiments have focused mostly on cropland and demonstrated the negative effects of topsoil reduction on crop production (Zhang F I G U R E 4 (a-e) Structural equation models (SEMs) for predicting productivity from several latent variables in different ecosystems.The blue solid lines indicate positive contributions from latent variables to productivity with coefficients provided next to the line.The red solid lines indicate negative contributions from latent variables to productivity with coefficients provided next to the line.Only the statistically significant coefficients are provided in the figure.The conceptual structure of the SEMs is shown in Figure S10.CFI, Comparative Fit Index; RMSEA, Root Mean Square Error of Approximation; TLI, Tucker Lewis Index.et al., 2021).Our results showed a clear positive association of topsoil depth with GPP in cropland, but such a relationship was not as strong as in grassland or shrubland.Grassland and shrubland root systems are often within 2-m depth compared to trees and thus more affected by topsoil properties.Moreover, grassland and shrubland were distributed mainly in drier regions (Figure S3), where topsoil depth was more strongly associated with ecosystem productivity.Under cropland and pasture, fertilization and irrigation may mask the effect of soil nutrients and water limitation due to deeper topsoils.It is noteworthy that soils under forest generally had shallow topsoils (mean thickness = 15 cm) but the highest GPP.There was no F I G U R E 5 Relationships between percent change in gross primary productivity (GPP) of a dry year (minimum precipitation) and a wet year (maximum precipitation) from 1986 to 2021 and topsoil depth in two climatic regions (arid and humid) and five ecosystems.Dashed horizontal lines indicate no change of GPP in specific conditions (y = 0).Dashed vertical lines indicate the mean topsoil depth in each region, which were used as threshold values to separate shallow and deep topsoils in Figure S13.
small.The GPP difference between deep and shallow topsoils and the GPP-topsoil depth relationship of the paired samples were not statistically significant across ecosystems and climatic regions.The change in GPP was generally positively related to change in topsoil depth, but such a relationship was not statistically significant.The ecosystem type and climatic regions were more related to GPP of the paired samples.One drawback of the paired comparison was that the topsoil depth and GPP data had different sample support [i.e., the length, area, or volume associated with a measurement(Goovaerts, 2014)].The topsoil depth was measured on point-based sampling locations, while the GPP was obtained from raster images with 30-m spatial resolution.Soil thickness has a large local-scale variation, and therefore the point-based measurements may not be able to represent its averaged depth in the 30-m raster pixel and may increase the randomness of the comparison at local scales, which can partially explain the negative relationship of GPP and topsoil depth at some watersheds.Additionally, other factors (e.g., ecosystem types and climatic factors) were more related to GPP, which may affect the GPP-topsoil depth relationship of the paired samples in different watersheds.Although topsoil depth was positively associated with the GPP in cropland, grassland, shrubland, and arid regions (Figure4; Figure S11 productivity in various ecosystems and different climatic conditions and results showed that such relationship was stronger in drier regions and grassland and shrubland.Deeper topsoil was also associated with improved resilience of ecosystem productivity to climatic extremes in these regions.The results can improve our understanding of the control of topsoil depth to ecosystem functioning and lead to a better representation of the role of soil in earth system modeling and climate modeling.It also provides evidence for natural resources management under climate change.However, limitations exist in the analysis and results interpretation.First, our analysis was based on existing observational data from long-term soil surveys, and we did not conduct any controlled experiments or explicitly select sampling locations based on edaphic, ecological, and climatic factors.The results and interpretation cannot completely eliminate multiple effects from other soil and environmental factors and solely investigate the effects of topsoil depth.Second, the large short-scale variation of topsoil depth and mismatch of the dimensions between soil data (point-based measurements) and GPP data (30-m resolution) may increase the randomness of the results at local scales.Third, this study investigated only the one-way effects of topsoil depth on ecosystem productivity ing.Future work is needed for creating national and global maps of topsoil depth.In addition, soil carbon stabilization and enhancement are considered as an important nature-based climate solution.Other soil factors (e.g., soil depth, bulk density) which are essential parameters to calculate SOC stock and directly related to soil erosion and compaction have been under-studied.Similarly, the role of other soil physical, chemical, and biological properties on ecosystem functioning under climate change should be further investigated.
Summary statistics of the paired comparison of GPP in different watersheds of five ecosystems in Figure 3.
TA B L E 1

GPP ∼ Depth binary × Ecosystem + SOC + Clay + Precipitation + Minimum temperature + (1| Watershed)
Unstandardized coefficients of fixed and random effect variables in the linear mixed-effect models with paired data.

Depth numeric × Ecosystem + SOC + Clay + Precipitation + Minimum temperature + (1| Watershed)
Topsoil depth was used as a binary variable (deep and shallow) and a numeric variable, respectively, in two models and interacted with ecosystem.The deep topsoil, cropland, shallow × cropland, and depth × cropland were considered as references for categorical variables.