Climate Change Dominated Long‐Term Soil Carbon Losses of Inner Mongolian Grasslands

Soil organic carbon (SOC) is the most critical component of global carbon cycle in grassland ecosystems. There has been growing interest in understanding SOC dynamics and driving forces of grassland biomes at various temporal and spatial scales. Up to now, estimates of long‐term and large‐scale changes in SOC of grassland biomes have been mostly based on modeling approaches and manipulative experiments, rather than direct measurements. During 2007–2011, we repeated 141 soil profiles of the sampling in 1963–1964 (up to 1‐m depth) to quantify the long‐term changes of SOC storage in the major grassland types of Inner Mongolia in order to tease apart the relative contributions of climate change and grazing. We found that SOC decreased in all soil types, except in the eolian sandy soils, from 1963 to 2007, with an average reduction rate of 1.8 kg C m−2 (~22.9% or 0.52% year−1) in the grassland biome of Inner Mongolia. We quantitatively clustered the soils into four groups using principal component analysis (PCA) and detected clear spatial dependency of the changes on climate and grazing. The climate change was responsible for 15.3–34.9% of the total SOC variations, whereas grazing intensity accounted for <9.5% of the changes. Our findings indicated that climate change, rather than grazing, was the primary forcing for the changes in SOC of Inner Mongolia grasslands. We presume that other driving forces, such as changes in nongrazing‐resultant wind erosion and atmospheric nitrogen deposition, might have played a role albeit their effects need to be further examined.


Introduction
Grasslands cover one fifth of the Earth's land surface and play an important role in the global terrestrial carbon cycle, especially under the changing climate (Hungate et al., 1997;Sala et al., 1996). Soil organic carbon (SOC) constitutes a major carbon component of grassland ecosystems, which also accounts for >80% of the global total (Parton et al., 1993) or 90% of China's grassland carbon stock (Ni, 2002). Minor changes in SOC may alter or even tip the carbon balance of grassland ecosystems (Davidson & Janssens, 2006;Trumbore & Czimczik, 2008). Additionally, the spatiotemporal changes of SOC, as well as the underlying mechanisms and feedbacks to climate, are jointly affected by the changing climate (Maia et al., 2019;Schimel et al., 1994;Scurlock & Hall, 1998) and human activities Derner et al., 2018;Larreguy et al., 2014;McSherry & Ritchie, 2013)-a research frontier in global climate change (Booker et al., 2013). Accurate estimation of SOC in grassland soils thus far is mostly based on modeling and experimental approaches that carry extremely high uncertainties at regional scale and over long-term periods (Groisman et al., 2018). A recent meta-analysis reported that increases in SOC content free from land-use changes accounts 48.1% of the studies (L. . Here, an overall increasing trend in SOC has been attributed to elevated precipitation and manure input through grazing (Goidts et al., 2009;Meersmans et al., 2009). Other studies have reported a decreasing SOC trend for the grasslands worldwide in recent decades due to multiple factors such as drought, over grazing, and other disturbances (Bellamy et al., 2005;Chapman et al., 2013;Fantappiè et al., 2011;Mestdagh et al., 2009;Schipper et al., 2007;Wang et al., 2003;Xie et al., 2007). These contradictory results reflect the substantial differences in vegetation and environment of grassland biome at global scale.
Several underlying mechanisms have been postulated to explain SOC dynamics for different geographical regions and at different spatial scales (Bell et al., 2011;Hartley, 2014), with a diverse suite of field and analytical methods among the studies (De Bruijn et al., 2012;Fang et al., 2010). For example, Bellamy et al. (2005) found that climate change was the primary cause for the significant loss of topsoil SOC in grasslands in England andWales during 1978-2003, although a modeling study of the same region concluded that only 10-20% of the SOC loss was due to the warming trend (Smith et al., 2007). McGovern et al. (2013), meanwhile, did not find any significant changes in SOC over a range of soil types and parent materials in the United Kindom in the past 40 years. They suggested that the soils were relatively resistant to large-scale changes in environmental factors such as temperature and grazing. Among the several potential reasons for these controversies is the lack of long-term in situ measurements of SOC, as well as the potential driving mechanisms at broader spatial scales (e.g., regional or national scale) Pineiro et al., 2009;Schipper et al., 2007;Wang et al., 2018).
With a total area of 792,000 km 2 , the grasslands of Inner Mongolia form the majority of the Euro-Asian steppe. In recent decades the grasslands have experienced a spatially and temporally variable climate (Figures 1c and 1d) (Gutman et al., 2020;John et al., 2018;Lu et al., 2009), as well as various types/ intensities of land uses, such as land conversion and grazing (John et al., 2016. Over the past half century, annual mean temperatures increased by 0.8-2.3°C throughout the Mongolian Plateau, while annual rainfall amounts in the region have decreased since 1998 (Groisman et al., , 2018Li & Wang, 2012;Liu et al., 2014;Lu et al., 2009). Meanwhile, human activities have been escalating. For example, grazing intensity witnessed a twofold to threefold increase from 1949 to 2010 (Figure 1b

10.1029/2020GB006559
Global Biogeochemical Cycles (J. Chen et al., 2015;Chen et al., 2018), which has been suggested to be the reason for the ecological degradation of 65% of the region's native grasslands (Chen et al., 2013;Li, 1997;State Environmental Protection Administration, 2008). These changes have been the catalyst for an increasing number of scientific endeavors on the functional changes of grassland ecosystems (e.g., species, production, energy balance, and SOC), as well as the independent and interactive contributions of climatic changes and land use (J. Zhao et al., 2015).
Several studies have reported changes in SOC dynamics over the past two decades in grasslands of China, including Inner Mongolia. These results appear to be inconsistent and inconclusive. Currently, SOC is generally assumed to be a significant atmospheric C source due to heavy grazing (Xie et al., 2007;Zou et al., 2007) and climate warming Zhao et al., 2015) or a neutral  or significant C sink owing to climate change (Huang et al., 2010;Piao et al., 2009). The differences among these estimates may be largely due to the data sources and methods used for estimating SOC. Ideally, analysis of soils from repeated sampling over a long period would provide the most direct and reliable answer. In 2007-2011, we repeated the soil surveys of 1963-1964 by the Chinese Academy of Sciences (CAS) in order to quantify the long-term spatial changes of SOC in the major grassland types of Inner Mongolia. A particular focus was to quantify the contributions of climate and grazing to the changes in SOC.
We hypothesized that climate change, rather than anthropogenic impacts (e.g., grazing), is the primary factor responsible for long-term changes in SOC in the native grasslands in Inner Mongolia. This hypothesis is detailed with several premises: (1) grazing may play a stronger role than climate change in regulating "fast change" variables such as above-and below-ground biomass in topsoils but has much less impact on organic carbon in deep soils; (2) grazing has been practiced on the Mongolian plateau for thousands of years with livestock as a key ecosystem component rather than an external force, implying that the ecosystem may have established adaptive mechanisms to cope with changes in livestock grazing at decadal scale (e.g., the east Mongolian Plateau in late 1100s-the Genghis Khan's era-had livestock density that was similar to the current level; Cui, 2011); and (3) regional temperature in Inner Mongolia has increased at a much higher rate than the global average since 1955 Lu et al., 2009) but with a highly variable spatial distribution of precipitation and extreme events such as dzuds and heatwaves Li & Wang, 2012). Unlike the changes in livestock density, climatic changes that occurred during our study period did not follow the long-term climatic changes before the warming trend in the most recent decades.

Original Soil Surveys
A field survey of physicochemical properties and nutrient status of the soils in major vegetation types of Inner Mongolia was conducted during 1963-1964 (hereafter the 1963 campaign) by the CAS team. The surveyed sites were strategically allocated by vegetation and soil type, with 163 soil profiles created and sampled in grasslands of Inner Mongolia. For each soil profile, the team recorded or measured the sampling location (description based on topographic maps), vegetation type, land use, soil horizons (including name, thickness, and major characteristics), organic matter content, and texture. Soil organic matter contents were determined using the Walkley-Black official procedure (Walkley & Black, 1934). SOC was calculated with a coefficient of 0.58 in this study.

Relocation and Soil Resampling
The soil profiles from the 1963 campaign were relocated and resampled during 2007-2011 (hereafter the 2007 campaign). The sampling locations of the 1963 campaign were recorded on 1:50,000 or 1:100,000 topographic maps, with an accuracy of 1 mi. Pinpointing their current locations involved several steps, using descriptions of each location, records of landmarks, topography, land use, vegetation type, and soil type from the 1963 database. First, we consulted several researchers who participated the 1963 campaign to approximate the sampling positions according to the historical records. Next, we took three soil cores (1 m) in the sites to generate descriptions of their vertical properties (e.g., layering, color, and texture by layer) for preliminary comparisons with those from 1963. Finally, when the properties of the soil cores at a site matched those of 1963 cores, a soil pit was created for renewed the sampling. One particular challenge in relocating and resampling was the trade-off between taking a greater number of samples at a smaller number of locations, versus taking fewer samples from a larger number of locations across the region. Based on our available human resources 10.1029/2020GB006559

Global Biogeochemical Cycles
(e.g., labor cost), traveling time, and the overall goal, we favored spreading our sampling efforts across the region to cover more soil/vegetation types, different climates, and diverse land uses. In the end, we were able to locate 141 sites that matched the 163 profiles. Seventeen of the historical profiles could not be relocated because of incomplete records in the database. Five profiles had clearly different visual features (e.g., layering, color, and texture) and were therefore excluded from this study.
To assess the accuracy of our site locations, we selected 30 sites from the 1963 database and compared their locations with other independent references, such as secondary use in different maps, topographic maps used in the 1963 campaign, grassland maps, and newer topographic maps that contain information on elevation, terrain, residential areas, traffic lines, and administrative boundaries that were developed before the 1980s (scaled at 1:100,000-1:250,000). We overlaid different maps to obtain a range that was likely to include the sampling locations from the 1963 campaign. The average (STD) distance between the repeated sampling sites during the two periods was 10.35 (4.77) m, with a maximum distance of 20.71 m (see details in the Supporting Information S1).
Each soil profile was photographed prior to sampling. Descriptions of the profiles to up to 1.0-m depth, or to the parent material horizon, were developed. Soil bulk density of each layer was sampled using a 100-cm 3 soil sampler, with approximately 1 kg of soils collected for each layer for laboratory analysis. Soil samples were oven-dried at 105°C for calculating bulk density. Roots were manually removed from the samples, with the remaining soil air-dried prior to physicochemical analysis. Rock fragments of >2 mm were separated with a 2-mm sieve and weighed. SOC content was measured using the Walkley-Black method, and soil texture (i.e., sand, silt, and clay proportions) was determined using the sieve-pipette method (Beuselinck et al., 1998); both of these methods were also used in 1963 campaign.

Climate and Livestock Data
The long-term climate data, including mean annual precipitation (MAP), mean annual temperature (MAT), accumulative temperature ≥5°C (T ≥5°C ), and surface evaporation (Evap) of 890 stations throughout China, were calculated from the records of the National Meteorology Data Centre of China and interpolated to 1 km × 1 km grids using the ordinary Kriging within the ArcGIS 9.3 software package. The accuracy evaluation for interpolated climate data appeared satisfactory (see details in the Supporting Information S1). The slopes (β) of MAT, MAP, T ≥5°C , and Evap from 1954 through 2007 were used to reflect the long-term changes (i.e., βMAT, βMAP, βT ≥5°C , and βEvap, respectively); decadal mean values of 1954-1963 and 1998-2007 were used to represent the average climate condition of the two sampling periods.
Average livestock population size by county during 1954-2007 was collected from the annual rural social economic investigation and assembled by the Inner Mongolia Bureau of Statistics. We used livestock density (LSK d ) to represent the mean grazing intensity, which was calculated as the livestock population divided by the grassland area and converted to the Dry Sheep Equivalent (DSE)--a unit representing the amount of feed stock required by a 2-year-old 50-kg Merino sheep to maintain its weight. Grassland areas were derived from the first to fourth Grassland Resource Inventory of Inner Mongolia. The slope of LSK d (βLSK d ) was used to represent the long-term changing rate of grazing intensity.

Data Analysis
SOC density (SOC D , kg C m −2 ) up to 1-m depth was estimated as follows: where n is the number of diagnostic horizons through the profile, BD i is bulk density (g cm −3 ), SOC i is organic carbon content (g kg −1 ), T i is thickness (cm)of the ith soil layer, and R vi (%) is the volume percentage of the fraction >2 mm (rock fragments). R vi was calculated from the weight percentage of rock fragments R mi as follows: where BD rf is the specific gravity of the rock fraction (i.e., 2.1 g cm −3 ) based on the average bulky density of the gravels ( "skeleton soil proportion" with diameters of 1-3, 3-5, and >5 mm. The rock fragment >2 mm was calculated from the skeleton soil proportion. Our preliminary examination of the relationship between particle sizes and skeleton soil proportion showed a clear exponential relationship. Consequently, we fitted the rock fragment content of >2 mm for each profile using the "Growth" function in EXCEL where y is the proportion of soil particles with a diameter greater than x and a and b are the empirical coefficients. For profiles that had less than two diameter grades, we used the following equation to estimate the rock fragment of >2 mm Soil bulk density (BD, g cm −3 ) was not measured in the 1963 campaign. We developed an empirical model between soil bulk density and SOC (%) (Bellamy et al., 2005;Yang et al., 2010) using the data set from the 2007 campaign to estimate the bulk density for 1963 We quantified the uncertainties of the BD estimation through bootstrapping resampling. The estimated BD and its uncertainties at the 95% confidence level were 1.2143 ± 0.2609 g cm −3 averaged across the 141 profiles, with considerations of uncertainties associated with the model parameters and prediction errors (see detailed descriptions of uncertainty analysis in the Supporting Information S1).

Quantifying the Spatiotemporal Change of SOC
The Shapiro-Wilk test was performed to examine the normality of all data sets. We found that SOC data were not normally distributed. Consequently, it was necessary to present the geometric mean, standard deviation (SD), and standard error (SE) of the data sets.
We also used other descriptive statistics (e.g., minimum, maximum, median, arithmetic mean, and frequency tables) to explore the statistical properties of SOC. We then carried out a suite of paired t test to determine if SOC by soil type was significantly different between 1963 and 2007. The variation in SOC in both campaigns appeared spatially dependent. Therefore, semivariance analysis was performed to explore the spatial distributions of SOC. After comparing the fitness of different semivariogram models, the exponential model was selected as the best one for this study.
To identify the major drivers of SOC dynamics under changes in climate and grazing activity, the data from 113 profiles taken from vegetation remain native grassland in both sampling periods were further explored through multivariate analysis. First, we applied principal component analysis (PCA) to explore the major environmental factors that may be responsible for the SOC dynamics. The SOC in both campaigns followed approximately a lognormal distribution and was therefore transformed with log 10 . The variables in the PCA included sampling location (latitude and longitude), topsoil texture (sand, silt, and clay content), climate baseline (MAT, T ≥5°C , MAP, and Evap) and long-term change rate (βMAT, βT ≥5°C , βMAP, and βEvap), and LSK d and change trend (βLSK d ) from 1954 through 2007. Second, the sites were clustered into groups according to PCA scores, administrative unit, and grassland utilization. Because the change in SOC followed a normal distribution and appeared spatially stationary, the general linear model was performed on different PCA groups to explore the contribution of climatic/human factors on SOC loss. The Type III Sum of Squares was adopted for analysis of variance (ANOVA) to evaluate the overall effects (i.e., independent and interactive effect) of potential drivers on the temporal changes of SOC. All statistical analyses were performed using R software version 3.3.1 (R Core Team, 2016).

Spatiotemporal Changes of SOCD
SOC decreased in all soil types from 1963 to 2007, except the eolian sandy soils, and with an average reduction amount of 1.8 kg C m −2 (~22.9% loss over the study period, or 0.52% year −1 ) ( be dramatic in the meadow soils, chestnut soils, and brown caliche soils, which accounted for 38% of the total land area and 60% of the grasslands (61.9 million ha) in Inner Mongolia. The geometric mean changed by −8.8 (p < 0.05), −2.5 (p < 0.001), and −0.9 kg C m −2 (p < 0.05) in meadow soils, chestnut soils, and brown caliche soils, respectively. In 1963, the SOC varied from a maximum of 31.4 kg C m −2 to a minimum of 1.0 kg C m −2 and followed a lognormal distribution (Shapiro-Wilk test for log 10 SOC, p value = 0.32), with a geometric mean value of 7.8 kg C m −2 (SE ± 0.2) and median value of 7.5 kg C m −2 (Figures 2a1 and  2a2). The SOC in 2007 showed a similar frequency distribution (Shapiro-Wilk test for log 10 SOC, p value = 0.15), but with a lowered geometric mean (6.0 kg C m −2 ) and variation (SE ± 0.2 kg C m −2 ; max = 37.6 kg C m −2 ; min = 0.7 kg C m −2 ; Md = 5.6 kg C m −2 ) (Figures 2b1 and 2b2). Interestingly, the change in SOC (ΔSOC) followed a normal distribution, with an average of −2.0 (±0.4) kg C m −2 ; 51% of the profiles' soil carbon reduction exceeded 1.0 kg C m −2 (Figures 2c1 and 2c2).
Among the 121 sampling profiles in the native steppes in the 1963 campaign, eight were tilled before being resampled in 2007, and 4 of 20 profiles that were cultivated before 1963 were abandoned and restored to grasslands. Altogether, only 113 profiles were retained as native steppe in both campaigns, wherein SOC decreased from 7.8 ± 0.2 kg C m −2 in 1963 to 5.9 ± 0.2 kg C m −2 in 2007 ( Table 2).
The semivariograms of SOC in both campaigns were successfully converged by applying the exponential model (Figures 2a3 and 2b3), with range of 385.1 km, sill of 65.1, and nugget of 0 in 1963, and range of 679.2 km, sill of 69.8, and nugget of 0 in 2007. The small nugget:sill ratio indicated strong spatial autocorrelation of SOC in both campaigns. However, the semivariogram of SOC change over time (i.e., ΔSOC) was featured with a pure nugget effect (Figure 2c3), suggesting that the change of SOC over the study period did not show any autocorrelation.

Factors Controlling SOC Changes in Native Grassland
The PCA on the dynamics of orthogonal variables indicated that the first two principal components accounted for 69.7% of the total variation, with SOC change possessing a weak relationship with both principal components, a slightly positive relationship with climate change (MAT, T ≥5°C , MAP, and Evap) and a neutral relationship with livestock density change (LSK d ). However, the PCA conducted on static environmental variables presented a higher cumulative variance for the first two principal components (78.7%), with change in SOC showing a negative relationship with temperature (MAT, T ≥5°C ) and livestock density and a Note. The geometric means were calculated because the SOC data were not normally distributed; the arithmetic means also were included for comparison with previous studies. The paired t test was performed. * Significant at p < 0.05. ** Significant at p < 0.01.

10.1029/2020GB006559
Global Biogeochemical Cycles neutral relationship with MAP and Evap. A strong relationship between SOC change and geographic location was found in both PCAs. The PCA based on dynamic and static environmental variables clustered the sites into four and three groups, respectively (Figures 3a and 3b), with both groups showing strong spatial dependency. Consequently, we combined these two PCA results and sorted the sites into four groups for further analysis (Figure 3c).
G1 was located in the northern part of Inner Mongolia, where the annual mean temperature was below 0 in 1960s. SOC changed slightly from 14.6 to 13.2 kg C m −2 (p > 0.05). However, significant declines were found for G2, G3, and G4, which were located in the central and southern regions of Inner Mongolia. SOC changed by −2.4 (p < 0.005), −1.8 (p < 0.005), and −1.9 kg C m −2 (p < 0.001), respectively (Figure 4a). Further analysis indicated that SOC loss increased with higher mean SOC in PCA group G2, G3, and G4. However, a negative relationship between ΔSOC and mean SOC was found in the PCA group G1 (Figure 4b).

10.1029/2020GB006559
Global Biogeochemical Cycles change was inversely correlated with βMAT in PCA groups G1 and G3, with the partial correlation coefficients of −0.505 (p < 0.05) and −0.450 (p < 0.05), respectively. The change in SOC also was significantly correlated with βMAP for G1, G3, and G4, with a partial correlation coefficient of 0.561 (p < 0.01), 0.534 (p < 0.01), and 0.344 (p < 0.05), respectively. The changes in SOC showed no significant relationship with the grazing intensity change rate for all four PCA groups. However, the fraction of variance in SOC change with the change in livestock density (βLSK d ) was 9.5% (p < 0.05) for G4 and < 3.0% for G1, G2, and G3.

SOC Dynamics Over the Past Four Decades
We resampled soil surveys of~45 years ago in major grassland types for direct and reliable estimates of the spatiotemporal changes of SOC in the grasslands of Inner Mongolia in this study. Our findings demonstrate that the SOC of grasslands in Inner Mongolia had declined significantly from 7.8 to 6.0 kg C m −2 during 1963-2007, averaging 0.52% year −1 . This loss is equivalent to 2.13 Pg C of SOC in 118.3-km 2 grasslands in Inner Mongolia, which amounts to one seventh of the cumulative emission of C produced by fossil fuel use into the atmosphere between 1980 and 2000 in China (Marland et al., 2013). Our results are comparable with previously reported values for Inner Mongolia. Wang et al. (2003) investigated

10.1029/2020GB006559
Global Biogeochemical Cycles the change of SOC stock between the 1960s and 1980s using the data from 236 and 2,473 soil profiles, respectively (i.e., the database from the First and Second National Soil Surveys of China). They found that SOC density in the grassland region of Inner Mongolia decreased from 2.15-5 to 0-2.15 kg m −2 , which equals a change rate of approximately 1% year −1, doubles our results. Xie et al. (2007) estimated the SOC stock in China and its change during the period from the 1980s to the 2000s using the data from the Second National Soil Survey of China (released in 1996 and2006). They reported a SOC reduction of 0.71 kg C m −2 in the grasslands of Inner Mongolia, which accounted for 9.5% of the cumulative soil C loss during 1980-2000, with a loss rate of 0.48% year −1 . In contrast, using data from 327 soil profiles from a soil survey and 275 soil profiles obtained from the Second National Soil Survey, Yang et al. (2010) described a relatively stable SOCD (6.63 kg C m −2 ) for the grasslands in northern China, including Inner Mongolia, with a change of 0.08 kg C m −2 from the 1980s to the 2000s. Based on soil carbon inventories extrapolated by satellite greenness measurements, ecosystem models, and atmospheric inversions, Piao et al. (2009) concluded that grassland soils in China generally acted as a significant C sink during 1980s-1990s.
There may exist several reasons for aforementioned discrepancies among the studies. First, data sources were different. Soil data from the First (1960s) and Second (1980s) National Soil Surveys were used in previous studies. These databases included few soil profiles from grassland regions: <300 soil profiles for the entire northern grassland region in the 1980s and even fewer in the 1960s. This small sample size prevents us from reaching a statistically robust conclusion. Second, these approaches may carry substantial uncertainty than direct measurements. Previous studies extrapolated the site-level records to regional values through spatial interpolations, development of empirical relationships with remote sensing reflectance or modeling. Third, the historical surveys always did not include bulk density and rock content. Empirical regression models based on current data were often used to estimate these key variables in calculating SOC, which would introduce additional errors. In this study, the bulky density (BD) in 1963 was obtained using an exponential model which may result in additional uncertainties. To assess its importance, we performed sensitivity analysis on the impacts of BD uncertainties on the SOCD change. We found that SOCD from the 2007 campaign (geometric mean = 6.0 kg C m −2 ) was significantly lower (p < 0.05) than the simulated BD (geometric mean: 6.6-9.3 kg C m −2 ) at different percentiles level (i.e., 5%, 10%…, 90%, and 95%) (Table S2). Our sensitivity analysis further confirms the differences in SOCD between 1963 and 2007. Cautiously, the estimation of BD can introduce some uncertainties; we should not rule out the possibility that modeled BD may results in very different SOC, especially under extreme conditions. Finally, previous studies were conducted for much larger regions in northern China that cover more types of grasslands, climate, and soils. Here, we used our resources to increase the sample size by focusing on three relative homogeneous grassland types in Inner Mongolia for a longer period (45 years). The average sampling size per land  Figure 3c). unit area was much higher than those in previous studies. With our results and other reported findings, it appears that SOC in the grasslands of Inner Mongolia had a sharp decrease from the 1960s to the 1980s, but this decreasing trend was much lower from the 1980s to the 2000s.

Driving Factors of SOC Change
Soil texture plays a key role in determining the long-term SOC dynamics in agricultural ecosystems (Li et al., 1994;Meersmans et al., 2009;Mulder et al., 2015). In this study, we found that the SOC loss was larger in soils with higher sand content. Specifically, we found that the changes in SOC in Inner Mongolia grasslands were positively correlated with clay content and negatively correlated with sand content (Figure 5). This supports the assumption of Wiesmeier et al. (2015), who concluded that carbon storage capacity of semiarid steppe soils depended highly on the silt and clay particles present in the soil. The fine component of soil texture tends to protect organic matter more effectively from decomposition than sandy soil, which is known to improve the passive organic carbon pool (Burke et al., 1989;Parton et al., 1993), while sandy soils are more prone to leaching losses.
Previous studies in Europe showed that SOC loss may also depend on the mean SOC content (Bellamy et al., 2005;Goidts et al., 2009). However, Chamberlain et al. (2010) found that the correlations between change in SOC and SOC content appeared spurious due to the regression to the mean effect. Kirk and Bellamy (2010) proposed a conceptual model for quantifying the changes in rates of SOC decomposition in relation to inputs from vegetation. They reported a negative relationship between the rate of change in SOC and mean SOC content in England and Wales. They further denoted that carbon-poor soil tended to gain carbon, whereas carbon-rich soil lost carbon, due to the facts that carbon-rich soils contain more stable carbon than carbon-poor soils (Schulze & Freibauer, 2005). Our current understanding of SOC dynamics is that SOC storage in an ecosystem is determined by two processes: C gains induced by plant litter (including litterfalls and roots) and C losses due to microbial decomposition. These two processes have different changing rates and could vary from place to place, resulting in a positive or negative change. In our study, we also found that the relationships between ΔSOC and mean SOC in PCA groups varied by landscapes across Inner Mongolia (Figure 4b). A positive relationship between the change in SOC and mean SOC content was found in the higher latitude region of Inner Mongolia (PCA G1), which was inconsistent with the results in PCA groups G2, G3, and G4. This suggests that the increase in temperature may raise the MAT above the biological threshold for plant growth in the higher latitude region. Vegetation inputs, in turn, could offset carbon loss caused by other factors (Ding et al., 2019;Giardina et al., 2014), particularly in high latitudes characterized by a wet and cold climate.
The relative contributions of climatic change and land use change to changes in ecosystem functions (e.g., ΔSOC) are among the top research priorities within the scientific community (John et al., 2016). Contradictory reports on the roles of climate and land use on SOC dynamics have been reported for England and Wales (Bellamy et al., 2005;Chamberlain et al., 2010;Smith et al., 2007), without clear-cut conclusions being drawn on this issue. In this study, we observed an overall increase in MAT by~1.6°C since the 1960s, an elevated spatial variation of MAP during the study period (also see Lu et al., 2009;Chen et al., 2013;John et al., 2016) and a synchronized decline of SOC ( Figure 6). By comparing the spatial patterns between ΔSOC and historical climatic changes, we accepted our alternative hypothesis that the rapid rise in MAT was the most significant variable responsible for the balance between the C gain and C loss. Although the relationship between ΔSOC and the climatic variables in our PCA analysis appeared weak, we found significant responses of SOC to climate warming and rainfall fluctuation when independently examined by PCA groups (Table 3). Clearly, our hypothesis that climate change dominates SOC dynamics needs to be tested spatially (i.e., location specific). Nevertheless, these findings highlight the remaining difficulties in understanding the response of SOC to climate change on a large scale, which can be synopsized in three aspects.
1. The overall mean change of SOC is not entirely representative of SOC dynamics because there exists a large spatial variation in climate. 2. Carbon losses in soils due to climatic warming could be offset by increases in net production input to the soil in the long term (Thornley & Cannell, 2001). 3. Our 45-year observation period may not be long enough to reflect an accumulation of climate-driven change in SOC (Goidts et al., 2009). Continued efforts are needed to accumulate even longer term field data to accurately test our hypothesis.

Global Biogeochemical Cycles
In the present study, only the long-term changes of MAT and MAP were included to represent the influences of climate change on the SOC dynamics for Inner Mongolian grasslands. Several relevant studies suggested that if long-term changes of an individual/particular season in temperature and precipitation were considered (e.g., spring-summer and autumn-winter periods), more in-depth mechanisms could be unraveled (Garcia-Pausas et al., 2007;Kühnel et al., 2019). For example, Kühnel et al. (2019) showed that increased autumn precipitations had led to decreases in SOC stock, whereas increased spring and summer precipitation caused increases in SOC stock in grassland soils of German. Other studies concluded that both spatial and temporal changes in net primary production and soil respiration were more dependent on the intra-annual variations of precipitation than on the annual amount of precipitation (Burke et al., 1997;Raich et al., 2002). Unfortunately, we are not able to confirm their findings due to limited data sources.
Grazing and other land uses have also been proposed as significant factors affecting SOC dynamics in many terrestrial ecosystems including drylands (e.g., Chen et al., 2013;Qi et al., 2017;Xun et al., 2018), with inconsistent and often contradictory results (Derner et al., 2018;Frank et al., 1995;Goidts et al., 2009;He et al., 2011;Reeder et al., 2004). In this study, livestock density (LSK d ) increased dramatically from 1963 to 2007 (Figure 7a). Grazed grasslands, which may have experienced grazing for +1,000 years, showed a significantly lower SOC (8.8 ± 1.0 and 6.7 ± 0.7 kg C m −2 ) compared to the ungrazed hayfields (13.7 ± 2.6 and 12.5 ± 2.4 kg C m −2 ) in both 1963 and 2007. This finding agrees well with Wang et al. (2003) and Xie et al. (2007), who attributed the carbon loss in Inner Mongolia to grazing. However, our multivariate analysis for different types showed no clear relationship between LSK d and ΔSOC. This implies that the historical management of grasslands may be more important than the current grazing pressure. In other words, grazing effects may be subject to long time delays and contribute little to simultaneous SOC dynamics (Conant et al., 2001). Additionally, it is possible that grazing may strongly influence aboveground vegetation and topsoil (e.g. roots and root exudation), while climatic change is more responsible for the deep, slow cycling of soil carbon .
Several studies have pointed out that soil erosion and atmospheric dry deposition may affect the thickness of the soil profile and thus alter SOC stock in arid and semiarid terrestrial ecosystems of north China (Yan et al., 2005;Zhang et al., 2014). It has been hypothesized, for example, that erosion could reduce SOC storage due to topsoil losses, and loess deposition could increase SOC storage at the deposition sites. However, in a forest-derived farm watershed, Liu et al. (2003) found that erosion reduced carbon emissions via dynamically replacing topsoil with subsoil that had lower SOC contents and higher passive SOC fractions, whereas it enhanced carbon uptake at the eroded sites by continuously taking away a fraction of SOC that can be replenished with enhanced plant residue input. As a result, soil erosion and deposition reduced CO 2 emissions from the soil into the atmosphere by exposing low carbon-bearing soil at eroded sites and by burying SOC at depositional sites. For our sites, soil erosion and dry deposition were virtually treated as dependent variables of grazing intensity. Consequently, their effects were partially retained or concealed for the overall effects of grazing intensity and climate change. We speculate herein that soil erosion may have posed less significant impacts on the dynamics of SOC stock, mostly because our study region experienced light soil erosion (i.e., 0.2 −1 kg m −2 year −1 ) (Yan et al., 2005). Here, the dry deposition rate was rather comparable to the erosion loss rate in most part of the Inner Mongolia steppe region, except lower reaches in southeastern Inner Mongolia where the sandland ranges predominate the landscapes (Li et al., 2000). Wang et al. (2000) reported an annual dry deposition of 35.2 kg m −2 year −1 as of 1994 that was resulted

10.1029/2020GB006559
Global Biogeochemical Cycles from an annual increment of 2.7-10% during the previous two decades in a southern Inner Mongolian steppe region. This explains in part the increases in SOC of the eolian soils detected in this study. In the resampling surveys, we did not find detectable changes in the layering and thickness from the soil profiles in most cases except some sites that had been seriously eroded by rainfall-resultant runoff. These sites were sparsely interspersed on the land surface. In other words, wind-induced erosion does not appear significant to lead to changes in SOC for the study region.

Conclusions
We demonstrated a cumulative reduction of approximately one-quarter of the total SOC stock, from 1963 to 2007, in the grassland biome of Inner Mongolia. We detected clear spatiotemporal dependency of the changes on climate and grazing, with climate change being the primary forcing for the variance in SOC and grazing being less significant. Specifically, we observed an overall increase in MAT, an elevated spatial variability in MAP, and a synchronized decline of SOC since the 1960s, implying that increases in temperature and the variable precipitation assumed major importance on SOC. Our findings suggest that the historical management of grasslands may be more important than the current grazing pressure; i.e., grazing effects may be subject to long time delays and contribute minor in SOC dynamics. Nevertheless, we highlight the remaining challenges in understanding the response of SOC to climate change on a large scale, albeit by recognizing that (1) the overall mean change of SOC is not entirely representative of SOC dynamics because of large spatial variations in climate and (2) the observation period mostly may not be long enough to reflect an accumulation of climate-driven change in SOC.