Geological setting and sampling

The accumulation of eolian material on the CLP initiated during the early Miocene22. In the Quaternary loess–paleosol sequence, paleosol units formed during interglacial episodes were characterized by slower eolian dust deposition, strengthened summer monsoonal rainfall, and enhanced pedogenesis relative to the glacials23. Large seasonal changes in precipitation and evaporation contribute to the wide occurrence of pedogenic carbonates across the CLP. Unlike the more humid southeastern CLP (e.g. typical sites including Lantian, Baoji), where extensive leaching completely decalcified the paleosols, or the more arid northwestern CLP (e.g. Huanxian, Jingyuan), where large amount of detrital carbonates remained due to insufficient rainfall24, our sampling location, the Luochuan section (35.76°N, 109.42°E) located in the central CLP (Supplementary Fig. 1), is characterized by the appropriate level of pedogenesis. The detrital carbonates are mostly absent in paleosol samples, whereas authigenic carbonates are still preserved (Supplementary Fig. 2). Before sampling, we trenched the soil profile (>1 m deep) to avoid contamination of regolith, and monitored the MS values of bulk paleosols, which were compared to published data25. Each sample was collected within a 10 cm interval. Paleosol samples used in this study mainly occur within the soil Bt/Bw horizons, characterized by subangular blocky structure with clay coatings attached on structural faces, and a 5–7.5 YR Munsell color.

Finely disseminated carbonates

In the field, the carbonates in our targeted paleosol samples have identifiable pedofeatures including finely disseminated cements filling the soil matrix, and carbonate hypocoatings (i.e. pseudomycelia) lining pores and tubules, which can be categorized as Stage 1 pedogenic carbonate26. Therefore, we term this carbonate fraction as FDC. SEM imaging further reveals that our sampled FDC are mainly composed of 5–10 μm-long needle fiber calcite (NFC), and clustered nanoscale calcite nanofibers (Fig. 1a–d). The SEM results are consistent with previous findings that NFC is commonly associated with smaller carbonate fibers such as nanofibers27, indicating similar formation conditions and a pedogenic origin. Importantly, NFC mostly appears between the humic horizon (O horizon) and subsurface C zone and vanishes with depth28. The location of NFC matches our field observations, which indicate the occurrence of FDC within the relatively shallow soil Bt/Bw horizon in comparison to the nodule-bearing Bk horizon. The origin of NFC is tied to both organic (fungal biomineralization) and inorganic activities (physicochemical precipitation)29,30. Whichever the case, it appears that the stable isotopic compositions (δ13C and δ18O) of NFC are indistinguishable from coeval calcite cements (rhombohedral calcite crystals)31, making bulk measurements of FDC well representative of all encapsulated calcite phases.

Fig. 1 Micromorphology and geochemistry of pedogenic carbonates. a Needle fiber calcites (NFC) in the inner channel of hypocoatings. b and c Randomly organized mesh of calcite nanofibers. d Sparsely distributed nanofibres lying on the surfaces of NFC. e Cross plots of Mn/Ca and Mg/Ca ratios of finely disseminated carbonate (FDC) samples in this study (red dots), compared with detrital carbonates from potential source regions (blue squares)32. Error bars represent standard deviations from multiple samples (n > 3). f Cross plots of Sr/Ca and Mg/Ca ratios of FDC samples in this study, compared with those of microcodium (green squares) in Holocene soils on the Chinese Loess Plateau33 Full size image

Incomplete dissolution of the detrital carbonate components (calcite and dolomite) from parent material of the loess (e.g. eolian dust) would bias the δ13C of the bulk soil carbonates. However, they appear to only impose very minor influence on our samples because we selected samples without the presence of any dolomite for subsequent analyses, as measured by the Fourier transform infrared spectrophotometry that is highly sensitive to carbonates24. Dolomite has a slower dissolution rate than calcite, and the final disappearance of detrital dolomite ensures the complete dissolution of detrital carbonates24. Further, we also measured trace elements of the FDC fractions in our paleosol samples. Previous work has shown that the concentrations of trace elements (e.g. Mn, Mg, and Sr) in pedogenic carbonates from the CLP are considerably lower than those of detrital carbonates derived from marine carbonate strata in source regions32, and are highly correlated with rainfall intensity33. In this study, we identified that the Sr/Ca and Mg/Ca ratios of the FDC fractions range from 0.72 to 1.48 mmol/mol (mean = 0.98 ± 0.31 mmol/mol) and 36–140 mmol/mol (mean = 81 ± 30 mmol/mol), respectively, both of which are generally consistent with those of microcodium (typical authigenic carbonate in paleosols) in Holocene soils on the CLP (Sr/Ca = 1.08 ± 0.68 mmol/mol, Mg/Ca = 48 ± 29 mmol/mol) (Fig. 1f)33. The Mn/Ca and Mg/Ca ratios of the FDC fractions are significantly lower than the carbonate samples from potential source regions of the loess— the deserts in Northern China (Fig. 1e)32, providing additional evidence for the pedogenic origin of our bulk paleosols carbonates.

Diagenetic alterations (i.e. the dissolution and reprecipitation of fine-grained carbonates) could also affect δ13C and other geochemical signatures. However, we expect this process to have minimal influence on our carbonate samples for two reasons. Firstly, the long-term trend of δ13C derived from the FDC fractions of our paleosol samples are almost identical with that determined from δ13C of coeval calcite nodules throughout the Pleistocene, although the nodule record has a much lower resolution (Supplementary Fig. 3). Calcite nodules are usually considered to be buffered against diagenesis11,18. Notably, the Pliocene paleosol formations in the CLP also show indistinguishable δ13C and δ18O values from mm-scale micrites to massive carbonate horizons34, suggesting that the FDC could potentially be used to reconstruct pCO 2 beyond the Pleistocene. Secondly, numerous studies of paleosol carbonates from the CLP suggest that their geochemical signals are pristine and record near-surface conditions. For instance, trace metal compositions of fine-grained pedogenic carbonates from loess–paleosol sections in the CLP appear to preserve local precipitation signals33 and show high-frequency variations on orbital timescales35.

Based on field observations, micromorphology, and geochemical characteristics, we conclude that the FDC fractions used in this study are of a pedogenic origin with minimal diagenetic alteration, thus suitable for paleosol-based pCO 2 reconstructions.

Resolving S(z)

Rearranging Eq. (1) and we obtain

$$S\left( z \right) = {{p{\mathrm{CO}}_2}} \times \frac{{{\mathrm{{\delta}}} ^{{\mathrm{13}}}{\mathrm{{C}}}_{\mathrm{{a}}} - \delta ^{{\mathrm{13}}}{\mathrm{{C}}}_{\mathrm{{s}}}}}{{{\mathrm{{\delta}}} ^{{\mathrm{13}}}{\mathrm{{C}}}_{\mathrm{{s}}} - {\mathrm{1}}{\mathrm{.0044}}{\mathrm{{\delta}}} ^{{\mathrm{13}}}{\mathrm{{C}}}_{\mathrm{{r}}} - {\mathrm{4}}{\mathrm{.4}}}}$$ (2)

which provides a mathematical solution for S(z). Carbonate and SOM were measured for stable C isotopes to obtain δ13C s , and δ13C r (see the “Methods” section). δ13C a was obtained from marine carbonate δ13C10. Atmospheric pCO 2 over the last 800 ky were derived from the ice core pCO 2 record36. Together they were used to compute S(z) of the Luochuan section over the last 800 ky. To estimate the uncertainties of the back-calculated S(z), we adopted Monte Carlo random sampling simulations to propagate errors from all input parameters (Supplementary Note 1). In brief, for each sample, we randomly sampled each parameter in Eq. (2) within its error range for 10,000 times and calculated a corresponding S(z) population. We report S(z) as the median values of the results from the Monte Carlo simulations, and define the associated uncertainties using the 16th and 84th percentiles of the simulated S(z) population (corresponding to the median value and the ±1σ error range in a standard Gaussian distribution).

When pedogenic carbonate precipitates deep in the soil profile where soil-respired CO 2 dominates over atmospheric CO 2 , back-calculated S(z) are inherently associated with sizable errors7. We therefore adopted the concept of an R factor, defined as the ratio between [CO 2 ] atm and S(z) (R = (δ13C s − 1.0044 × δ13C r − 4.4)/(δ13C a − δ13C s )), as a screening criterion7. When atmospheric CO 2 constitute a minor portion of soil CO 2 (R < 0.3 according to ref. 7), the samples are thought to be dominated by soil-respired CO 2 thus discarded for subsequent analyses. Twenty-two paleosol samples passed this criterion, representing 78% of all the samples analyzed for S(z). The calculated median S(z) levels over the last 800 ky show variations from 356 to 815 ppm, with uncertainties on average of +96/−75 ppm (Fig. 2a).

Fig. 2 The magnetic susceptibility (MS)–S(z) correlation of the Chinese Loess Plateau. a Back-calculated S(z) values plotted against MS from Luochuan paleosol samples spanning the last 800 ky. Error bars represent 1σ errors resulted from Monte Carlo simulations (Supplementary Note 1). The best fitting line, as well as 95% confidence interval are indicated by red solid and dotted lines, respectively. b Histogram of the mean relative differences (χ) between calculated and ice core pCO 2 , when the sample numbers in the training group (n) is set to 10 Full size image

MS as a proxy for S(z)

Soil-respired CO 2 comes from two major sources—the respiration of autotrophs, such as plant roots, and heterotrophs, such as the soil microbes mediating the oxidation of SOM. Studies of modern soils have shown that S(z) as a measure of soil productivity, is strongly correlated to climate parameters, such as temperature, precipitation, and evapotranspiration37. Notably, among carbonate-bearing soils in semi-arid to arid regions, such as the CLP, S(z) is mainly controlled by precipitation18,37. To quantitatively constrain S(z), here we explore a potential proxy sensitive to precipitation changes—soil magnetic susceptibility (MS).

The MS of bulk soil samples is predominantly controlled by the presence of ultrafine magnetite38,39, which is efficiently produced in situ under well-drained soils with alternative wet (Fe2+ produced) and dry (Fe2+/Fe3+ oxide precipitated) cycles. The ultrafine ferromagnetic grains are either directly produced by magnetotactic bacteria, or formed by inorganic precipitation, mediated by the production of Fe2+ through iron-reducing bacteria40. Similarly, modest precipitation increase would enhance microbial activities, the formation of magnetic iron oxides and consequently soil MS, but MS decreases under excessive rainfall (MAP > 1000 mm) as the production of ferrimagnet ceases41. MAP > 1000 mm probably never occurred in the Pleistocene history of the central to northern part of the CLP42. As a result, MS is regarded as a sensitive paleoclimate proxy of pedogenesis intensity and paleorainfall40, and has been widely used to address the evolution history of the East Asian summer monsoon23,43.

Intensified precipitation would facilitate the formation of magnetic iron oxides (i.e. MS), and increase soil productivity and S(z). Because rainfall regulates soil MS and S(z) in a similar manner on the CLP, these two parameters are expected to be correlated. Indeed, when the two parameters from the Luochuan samples covering the last 800 ky were plotted against each other, we observe a statistically significant correlation (Fig. 2a):

$$S(z) = 2.66( \!\pm \!0.44) \times {\mathrm{{MS}}} + 114.9( \!\pm \!71.1)(r^2 = 0.64,{p} \, < \, 0.0001)$$ (3)

The fitting is done using the least-squares linear regression model, and the uncertainties on the slope and intercept represent 1σ standard errors. Extrapolating this relationship to the entire Pleistocene will enable sample-specific S(z) using measured MS, and consequently the calculation of pCO 2 from paleosols (MS–S(z) approach). By providing new constraints on S(z), this MS–S(z) approach potentially improves the paleosol method for reconstructing pCO 2 .

To validate the MS–S(z) approach, we next recalculated the S(z) of our 800-ky paleosol samples using the MS–S(z) approach, and compare the derived pCO 2 to those documented in ice cores. Because the determination of the MS–S(z) relationship relied on the ice core CO 2 record, to avoid circular reasoning, we adopted a resampling technique. Specifically, the 22 samples were divided into two subsets: a training sample group and a test sample group. A MS–S(z) regression was established from the training group, and used to estimate S(z) for the test group, the resulting pCO 2 of which were then compared to the ice core data. For a given training sample number n, a bootstrap sampling of 1,000,000 times was performed (each time, we obtained 22 − n pCO 2 values). We also varied the sample numbers in the training group from n = 10 to n = 21. During each run, the mean relative difference χ between the calculated pCO 2 and those from the ice core data was calculated:

$$\chi = \frac{1}{{\left( {{\mathrm{22}} - n} \right)}} \times \mathop {\sum}\limits_{i = 1}^{22 - n} {\left( {1 - \frac{{{\mathrm{{Calculated}}}\;p{\mathrm{{CO}}}_2}}{{{\mathrm{{Ice}}}\;{\mathrm{{core}}}\;p{\mathrm{{CO}}}_2}}} \right)}$$ (4)

The statistical distributions of χ (Fig. 2b and Supplementary Fig. 4) show that, across different training sample number n, χ is clustered around 0 with >70% of the data points falling within ± 10% difference. This consistent pattern among different n values strongly supports our MS–S(z) approach for pCO 2 reconstruction over the last 800 ky.

Reconstructing early Pleistocene pCO 2

The refined soil carbonate-pCO 2 approach was also applied to the samples from the lower section of the Luochuan profile spanning the early Pleistocene, beyond the available ice core records. Measured MS and stable carbon isotopes for those early Pleistocene samples were used to calculate S(z) and pCO 2 , assuming that the MS–S(z) relationship remains (Eq. (3)). We calculated the uncertainties on pCO 2 by propagating the uncertainties associated with the analyses of all input parameters and the MS–S(z) regression equation using PBUQ, a published MatLab-based program that adopts Monte Carlo random sampling simulations for uncertainty analysis7. We report the median values of the calculated pCO 2 , and define the uncertainties as the 16th and 84th percentiles (as ±1σ), and the 2.5th and 97.5th percentiles (as ±2σ) of pCO 2 distributions from Monte Carlo simulations (see the “Methods” section).

The calculated S(z) during ~2.6–0.8 Ma vary between 276 and 644 ppm, with errors (1σ) ranging from 76 to 113 ppm (Fig. 3c). The R values (pCO 2 /S(z)) are generally higher (>0.3) during the early Pleistocene, suggesting that these samples are suitable for pCO 2 reconstructions. This newly established interglacial pCO 2 records show variations from 183 to 292 ppm (averaged median pCO 2 levels for each interglacial episode) during 2.6–0.9 Ma (Fig. 3d). Except for some data points centered around the Pliocene–Pleistocene boundary (2.6–2.5 Ma) and the MPT showing relatively higher pCO 2 exceeding 300 ppm, our paleosol-CO 2 estimates document overall low early Pleistocene pCO 2 levels similar to those over the last 800 ky (Fig. 3d). The reconstructed pCO 2 values have an 1σ error range of +92/−82 ppm when averaging the errors for all individual results.

Fig. 3 Stable isotopes and calculated S(z) and pCO 2 . a δ13C values of soil organic matter (δ13C SOM ). b δ13C values of finely disseminated carbonates (δ13C c ). c S(z) estimates based on the magnetic susceptibility (MS) proxy. d Reconstructed early Pleistocene pCO 2 (median levels of pCO 2 distributions, 2.6–0.8 Ma) and ice core CO 2 record (0.8–0 Ma)36. Horizontal gray line shows the pre-industrial pCO 2 level (280 ppm). Errors associated with δ13C c and δ13C SOM are standard deviations of all measurements within the same paleosol unit (n > 3). Error bars related to pCO 2 represent the 16th and 84th percentiles (1σ) based on PBUQ Monte Carlo simulations7, whereas errors associated with S(z) were calculated based on Gaussian error propagations (Supplementary Note 1) Full size image

To evaluate the role pCO 2 played for the MPT, we performed a detailed comparison of MPT–pCO 2 with those before and after the MPT (i.e. pre-MPT and post-MPT), the durations of which were considered to be 400 ky equivalent to the MPT itself (i.e. 0.8–0.4 Ma for post-MPT and 1.6–1.2 Ma for pre-MPT). Our paleosol-based CO 2 reconstructions provide 16 estimates during the MPT and 21 pre-MPT. pCO 2 estimates during the MPT and pre-MPT were then converted into factor changes in CO 2 by normalizing pCO 2 to the mean post-MPT pCO 2 from ice cores, since our paleosol–CO 2 estimates are indistinguishable from ice core data (Fig. 2b and Supplementary Fig. 4). We calculated the probability density functions (PDFs) and cumulative distribution functions (CDFs) of the factor changes in CO 2 among different time periods (i.e. MPT/pre-MPT, MPT/post-MPT, pre-MPT/post-MPT) (Supplementary Note 2). As indicated by both absolute values and the factor changes, pCO 2 level during the MPT (mean = 269 ± 38 ppm, 1σ = +105/−94 ppm) was ~15% higher than that of the post-MPT (mean = 231 ± 37 ppm), and ~20% higher than that of the pre-MPT (mean = 217 ± 25 ppm, 1σ = +73/−67 ppm) (Fig. 4a). Results from the PDFs and CDFs also confirm this distribution (Fig. 4b, c).