For this study, we used the total soil and SOM fractions from Lugato et al. (2021). Briefly, 400 representative soil samples from the Land Use and Coverage Area frame Survey (LUCAS; Orgiazzi et al., 2018) across the European Union and United Kingdom were selected for SOM fractionation (Supplementary Figure S1). These were all mineral (soil organic C < 120gCkg⁻¹) topsoils (0–20cm). Size fractionation was conducted to separate POM (>53μm) from MAOM (<53μm) by wet sieving after mechanical dispersion of aggregates by shaking in sodium hexametaphosphate with beads, as described in Cotrufo et al. (2019). Briefly, for each sample 5g of 2mm sieved and oven dried soil was shaken in dilute (0.5%) sodium hexametaphosphate with beads for 18h to break the aggregates. The dispersed soil was then rinsed onto a 53µm sieve, with the fraction remaining on the sieve being collected as POM, and the fraction passing through (<53µm) collected as MAOM. This fractionation approach, while convenient for the high throughput, defines POM by size (>53μm), thus small amounts of very fine POM could be recovered in the MAOM fraction, and OM associated to coarse minerals (e.g., sand) is recovered as POM. A separation of POM by size and density, when feasible, is advisable. All bulk soils, MAOM and POM samples were pulverized, oven dried at 60°C and analyzed for C and N concentration in an elemental analyzer (LECO TruSpec CN). Inorganic C was removed by acid fumigation (Harris et al., 2001) before elemental analyses, in the few samples where it was present. Due to a variety of reasons (non-readable labels, broken sample containers during transportation and general quality checked after fractionation), 349 out of 400 were retained for MIR and statistical analysis.

For each soil sample we extracted data on particle size distribution (% clay, silt, and sand), pH, available phosphorus (p mgkg⁻¹), extractable potassium (K mgkg⁻¹) cation exchange capacity (CEC cmol (+) kg⁻¹), and land use types from the LUCAS database (http://esdac.jrc.ec.europa.eu/content/lucas-2009-topsoil-data). Bulk density (BD) was indirectly estimated using a pedotransfer functions (Hollis et al., 2012). Values and basic statistics of physicochemical soil properties are summarized in Table 1.

Long-term mean annual temperature (MAT, °C) and annual precipitation (p, mm) were obtained from the high-resolution WorldClim datasets (version 2), which has average monthly climate data for minimum, mean, and maximum temperature and for precipitation for 1970–2000 (http://worldclim.org) at 30s resolution. Using these data, the aridity index (AI) (UNEP, 1997) was calculated as the ratio of mean annual precipitation (MAP) to potential evapotranspiration (ET). We used the aridity index to separate our soils into different climatic zones: arid (0.03–0.2), semiarid (0.2–0.5), dry sub-humid (0.5–0.65), wet sub-humid (0.65–1) and humid (>1.0).

The 349 bulk soil samples were air-dried, ground to fine powder, and scanned undiluted (neat) in the mid-infrared spectral region (4,000–400cm⁻¹). The spectra were obtained using a Digilab FTS 7000 spectrometer (Agilent Technologies, Walnut Creek, CA) with a KBr beam splitter, and a Peltier-cooled DTGS detector. The samples were analyzed in diffuse reflectance, and each spectrum consisted of 64 co-added scans at 4cm⁻¹ resolution. Scans were conducted in duplicate for each sample, and subsequently, the average spectrum was calculated using GRAMS/AI software (Thermo Scientific).

Partial least Square (PLS) regressions were performed to estimate C and N pools, clay, pH, C/N ratio and soil texture (clay and sand content). Several mathematical pre-treatments were tested, including the Savitzky–Golay function, the second derivative, standard normal variate (SNV), and multiplicative scatter correction (MSC) to reduce baseline variations and/or enhance spectral features. The best improvement was reached using Savitzky–Golay functions with a first derivative using three point smoothing (Savitzky and Golay, 1964). A non-linear Kernel partial least squares (PLS) regression model was used for building the predictive models. During model calibration, a priori outlier detection techniques were adopted to remove influential outliers by plotting the residual X- and Y-variances against leverages and critical thresholds as calculated from Hotelling's T 2 and Q-residuals statistics. The fraction of outliers removed from the models varied from 2 to 2.5%. The PLS calibrations were validated using a random cross-validation approach. The cross-validation was performed by randomly dividing the calibration set into multiple segments using the default setting of the software (Westad and Marini, 2015). The spectral analysis and the development of PLS models were carried out in the Unscrambler chemometric software V.10.3 (CAMO, Norway).

The accuracies and the robustness of the PLS models were evaluated for calibration and prediction using the root mean square error (RMSE) and root mean square error of cross-validation (RMSEcv), determination coefficient (R ²) and the ratio of prediction to deviation (RPD), which was calculated as the standard deviation divided by RMSE. The RMSE is a measure of the difference between the predicted value and measured value, while the R ² indicates how well the model calibration and validation fits the measured data. RPD is a dimensionless parameter that represents the ratio of the standard deviation to the RMSE.

The model performance was classified according to the following RPD intervals (Viscarra Rossel et al., 2006; Bellon-Maurel et al., 2010): RPD ≥2 to indicate the model performance was very good to excellent and reliable for quantitative predictions; 1.8 < RPD <2 to indicate satisfactory performance for quantitative predictions, but may be further improved by other approaches; RPD <1.8 to indicates fair predictions and RPD <1.4 poor or insufficient for most applications.

The database was analyzed by Principal component analysis (PCA) and subsequent redundancy analysis (RDA) to elucidate the importance of explanatory variables such as land use types, climate and soil properties on soil spectral properties. Before executing the PCA, the averaged spectra for each soil were baseline corrected and mean-centered. In the RDA, the most significant discriminating variables on spectral bands were selected based on PCA loadings, which was carried out with CANOCO version 5.01 (Microcomputer Power, United States). A correlation matrix provided in the supplementary material was calculated to determine the strength of the relationships between soil properties using the corrplot package in R. A p-value of less than 0.05 was considered statistically significant. The Pearson’s correlation (r) analysis between MIR absorbance values and C and N content in SOM, MAOM and POM were represented using contour/heatmaps. Correction for baseline offsets was applied to the spectra, and the results (r) were linearly interpolated into a grid to smooth the contour/heatmaps (Sigmaplot v.11).

We performed PCA using the full MIR spectral range of the sample set to determine the influence of climate and land use on the spectral response of bulk soils. The first component explained 69% of the variance, whereas the second and third components accounted for 17 and 5%, respectively. In general, the sample score distribution shows that the variation in the spectral signatures of soils was associated with differences in the aridity index (AI) and land use types (Figure 1). In the PCA analysis, samples from arid climates were separated from soils under more humid conditions (Figure 1A). For land use types, only woodland soils had a noticeably different pattern from other land uses, which tended to overlap with each other (Figure 1B). Furthermore, the sample score distribution shows that the variation in the spectral features was also associated with differences in C content. Samples were separated from soils that exhibited low and high SOM C content, however, some samples with contrasting C content also exhibited similar spectral composition (Supplementary Figure S2).

To further explore these differences, we used the PCA loadings and score distribution to reduce the dimensionality of the spectral data, as well as to identify the spectral features contributing to the maximum variance within the sampling set. The loading plot (Figure 1C) showed that various spectral bands from organic compounds had positive loadings in PC1. The spectral differences were associated with higher absorbances in the 3,000–1,500cm⁻¹ region, attributed to the presence of C–H stretching vibrations and aromatic groups near 1,600–1,570cm⁻¹ and in the 2,900–2,940 cm⁻¹ region. A pronounced band characterized the second loading in the 1,200–1,042cm⁻¹ corresponding to C-O stretching of polysaccharides. It should be noted that the band of polysaccharides may overlap with mineral compound bands at 2000–1800cm⁻¹ and 1,030cm⁻¹, which can also be attributed to Si-O groups of quartz and clay minerals, respectively. Furthermore, we found positive loading bands at 2520^–1cm and 3,620cm⁻¹ corresponding to carbonates and stretching vibrations of OH groups in lattice clay minerals, respectively.

The relationship among soil properties, land use type and climate, along with a selected subset of spectral bands, were identified using redundancy analysis (RDA) (Figure 2). The bands in the RDA were selected in order of the highest loading per component (Figure 1A). Bare land and shrubland were excluded from this analysis due to the small number of samples within these categories. The results show that C accumulation in the MAOM fraction (ƒ MAOM) was particularly high under agricultural conditions with inherently low SOM C (Figure 2, Supplementary Table S1). Furthermore, clay was found to significantly increase the C storage in the MAOM fraction (Figure 2, Supplementary Figure S3). These soil features were associated with the absorption bands at 2,510cm⁻¹ and 3,624cm⁻¹ (Figure 2), corresponding to carbonate and clay minerals, respectively. In more humid conditions, woodland and grassland were positively correlated with SOM C, sand, and C and N content in the POM fraction. In addition, these conditions show a substantial increase in various functional groups attributed to aliphatic C (3,000–2,800cm⁻¹), carboxyl C=O (1700–1800cm⁻¹), aromatic C (1,600–1,500cm⁻¹) and amide N-H (1,550–1,500cm⁻¹).

In large and complex datasets, the exploration of spectral variations associated with local conditions might be used to identify uninformative MIR regions. In this way, redundancy and collinearity are reduced so as to increase the model accuracy. Smoothed heatmaps show that the Pearson’s correlation coefficient (r) between absorbances at MIR wavelengths and concentrations of C in SOM, MAOM and POM is land use-dependent (Figure 3, Supplementary Figure S4). Heatmaps illustrate that overall, the position and the magnitude of correlation between spectral bands and the C content in SOM, MAOM and POM differed among woodland, grassland and cropland (Figure 3). For example, POM fraction heatmaps showed a stronger correlation in the 2000–1,000cm⁻¹ region for grassland and woodland (Figure 3C). In contrast, this region was almost indiscernible for cropland. In addition, for woodland, evident differences were observed in the organic matter and mineral absorbances in MAOM fraction, especially in MIR region (2000–1,000cm⁻¹). Furthermore, correlation analysis also shows other spectral patterns that are related to soil organic or inorganic C pools (Supplementary Figure S4). We found high correlations between C in SOM, MAOM and POM, and aliphatic bands around 2,850–2,920cm⁻¹ for grassland and woodland, but less so for cropland. In addition, woodland exhibited a low correlation with the Si−O groups in quartz near 1800cm⁻¹ for SOM C and MAOM C (Supplementary Figure S4). For inorganic C, the band at 2,525–2,510cm⁻¹ had a strong correlation with CaCO₃ content in the grassland and cropland, but less so in the woodlands. The lattice clay absorbance bands between 3,650 and 3,595cm⁻¹ were positively correlated to clay content as expected, but slightly negatively correlated to organic C and POM C in soils from every land use type.

The performance ability of calibration models is shown in Table 2. Overall, low values of RMSE and high R ² indicate better predictive performance, however, RPD is used to standardize RMSE values against the SD, as RMSE is a dependent concentration parameter. For C pools, the best PLS-MIR model was developed within the 2000–400cm⁻¹ spectral range, which included the major absorption bands related to C pools previously identified in Figure 3. For C in total SOM and MAOM (Figure 4A; Figure 4B), good and reliable estimates were obtained using a global calibration (R ² = 0.80–0.79; RPD = 2.23–2.17, respectively). However, the cross-validation for C in SOM tended to underestimate the observed concentrations, thus a large number of data points were distributed below the 1:1 line showing a high dispersion at C content over 50gCkg⁻¹ soil (Figure 4A). In contrast, poor prediction accuracy was obtained for POM C showing a curvilinear trend (R ² = 0.61; RPD = 1.60). The model had a clear tendency to underestimate the prediction when the measured C content was higher than 20gCkg⁻¹ soil (Figure 4C).

For N pools, the use of the full MIR range demonstrated a better performance than models constructed using a wavelength selection (Table 2). Although the prediction of N in total SOM and MAOM was less satisfactory than for C, good results were also achieved for these variables (R ² = 0.64–0.51; RPD = 1.88–2.08, respectively) (Figures 5A,B). Regarding POM N (Figure 5C), a poor prediction was obtained and insufficient for most applications, showing similar behavior to POM-C (R ² = 0.47, RPD = 1.37).

The distribution of C and N in MAOM and POM fractions is highly controlled by predictors such as C/N ratio, pH, clay and sand (Figure 2). We compared the PLS model performance for the aforementioned soil properties (Table 2; Figure 6). Predictions for soil C/N (Figure 6A) were not improved (R ² = 0.67; RPD = 1.75) when compared with the individual predictions for C and N in SOM. For pH (Figure 6B), the best PLS-MIR model was developed using the 400–4,000cm⁻¹ spectral range, which showed a good prediction performance (R ² = 0.76; RPD = 2.06). Regarding the soil texture variables (Figures 6C,D), similar accuracies and excellent estimates were achieved for both clay (R ² = 0.81; RPD = 2.28) and sand (R ² = 0.83; RPD = 2.41).

Understanding and capturing the spatial variability of soil spectral signatures at regional to continental scales is relevant to developing a robust global calibration strategy and to subsequently improving the predictive ability of soil C models. Our results showed inherent soil spectral heterogeneity associated with climate and land use (Figure 1). Likewise, climate and land use are widely accepted as the most significant drivers of C content when we evaluate their distribution at large geographical scale (Rial et al., 2017; Ogle et al., 2019). Our data show that soils under drier conditions tend to have a relatively high soil pH (Figure 2), which is not surprising given the existence of calcareous materials characterized by the presence of carbonates at the 2,513cm⁻¹ band, concurrent with low organic matter content (Asgari et al., 2020). In semiarid environments, the decline in primary productivity and increase in litter decomposition due to higher temperatures and intermittent precipitation has resulted in less C storage when compared to colder and more humid conditions (Guo and Gifford, 2002). In fact, we observed low POM C content in agricultural systems under semiarid conditions (Figure 2). In addition, we observed a higher proportion of C in the MAOM fraction (f MAOM) and clay content. The latter responses are consistent with the correlation observed for the clay band at 3,620cm⁻¹, as well as the presence of organic matter absorbance around the 1,100cm⁻¹ band of polysaccharides (Rasmussen et al., 2018).

MIR spectroscopy has been used for qualitative analysis of various soil properties because it contains extensive information on the functional moieties of organic matter (Ramírez et al., 2020; Dudek et al., 2021). Additionally, it is a useful tool to explore spectral patterns associated with soil C stocks (Stumpe et al., 2011). For example, as in many previous studies, we also found aliphatic CH spectral bands at 2,926cm⁻¹ (Figure 2), which could be used to rapidly estimate soil C, and N concentrations of northern cold region soils (Matamala et al., 2017). Likewise, we also found distinct bands in the region between 1,510 and 1,230cm⁻¹, which have already been reported in association with the decomposition of forest litter in soils (Haberhauer et al., 1998). It should be noted, however, that band assignments are not absolute but rather are an approximation in the sense that a band for a particular compound is assigned at a specific wavenumber.

The MIR spectral bands can be used as valuable data to monitor major chemical changes from different land use systems. It is interesting to note that although MIR is not a high-resolution tool, the information provided by MIR spectroscopy was able to elucidate spectral patterns across different C pools. For example, the spectral patterns between POM and MAOM involved variation in absorbances between the 2,800 and 1,000cm⁻¹ region (Figure 3). This spectral region is characterized by silicates (Si-O vibrations) and different C functional groups, including aliphatic C–H, carbonyls, aromatic C and C-H stretching typically found in carbohydrates (Parikh et al., 2014). This analysis revealed that POM displays high spectral similarities to SOM C in woodland systems, while cropland and grassland did not. Furthermore, in the cropland systems, we found similar spectral patterns between SOM and MAOM, as well as a substantial decrease of organic C absorbances in the POM fraction. This is consistent with the relative distribution of C (f MAOM) between these fractions (Supplementary Table S1) with a higher amount of C in POM in woodland soils with overall higher C content and lower disturbance as compared to grassland and cropland soils (Lugato et al., 2021).

We achieved accurate global predictions using MIR coupled with PLS spectroscopy for C and N pools at a continental scale. In the absence of some of the wavelengths containing unrelated information, the wavelengths 2000–400cm⁻¹ had a substantial influence on the PLS model for the quantification of C in SOM and MAOM. As often reported for soil C, the multiple and small bands around 830–400cm⁻¹, as well as the fingerprint region (1,500cm⁻¹–650cm⁻¹) contribute to improve the predictive ability of PLS models (Cécillon et al., 2012; Knox et al., 2015). Our results showed that MIR produced accurate predictions to determine C and N in SOM and MAOM (Figure 4). Typically, the prediction accuracies reported for total soil N are similar to the accuracies of the SOM C models. This is mainly because most of the N is bound in the SOM, and thus, soil N is often tightly correlated to total soil C (Schirrmann et al., 2013). C and N in POM were less accurately predicted than other C or N pools (Figures 4C, 5C). For example, our prediction yielded a low accuracy for POM C (R ² = 0.61), indicating that the PLS calibration could not capture the variability over a wide range of 0.46–54.9gCkg⁻¹. Although previous studies have reported satisfactory accuracies for C prediction in the POM fraction using MIR spectroscopy (R ² = 0.71–0.83), they were carried out in the 0.2–16.8gCkg⁻¹ range of concentration (Janik et al., 2007; Baldock et al., 2018). However, recent studies using heterogeneous soil samples at a country level similarly reported that SOM C and MAOM C were better predicted than POM C (Sanderman et al., 2020). These results can be attributed to the fact that POM is more heterogeneous than MAOM when isolated by size, since it includes both a light fraction not associated with minerals, as well as heavy coarse organic matter associated with sand-size minerals (Lavallee et al., 2020). Therefore, the spectra are subject to a high degree of variability due to the various decomposition stages of plant residues, and to the range of contribution of compounds of microbial origin which adversely impact the calibration accuracy. Furthermore, the light component of POM mostly consists of partially decomposed plant material, so differences in the vegetation cover (woody plant cover vs grasses, for example) affect the chemical character of POM (Soucémarianadin et al., 2019).

We found that C estimations in SOM were more reliable than those provided by studies that used NIR spectroscopy in the same large-scale European dataset (Nocita et al., 2014; Rial et al., 2017; Ward et al., 2019). For example, the selection of sample subsets by k-Means clustering achieved less accurate predictions (R ² = 0.67; RPD = 1.74) (Ward et al., 2019). In particular, MIR spectra contain fundamental bands for several organic and mineral moieties, while soil vis-NIR ranges are largely non-specific, consisting of weak, broad and overlapping absorption bands (Reeves, 2010). Furthermore, we found that MIR was highly sensitive to soil mineral composition, and appeared well equipped to determine pH and soil texture, such as clay and sand using PLS models (Figure 6). It is important to considerer we used a size separation of POM and MAOM. Therefore, as mentioned above, the organic matter associated to sand-size particles was considered POM, and could have contributed proportionally more to this fraction in coarse-textured soils. However, despite this limitation, we suggest that the influence of sand content was probably low for the prediction of C and N in POM. MIR spectra present greater signatures in quartz and clay minerals, which provide more robust calibrations to predict clay and sand content (Viscarra Rossel et al., 2006; Reeves, 2010). Therefore, in the presence of quartz-rich sand content, we clearly would expect more accurate predictions for POM than those achieved in this study. Another limitation of our study was the negative predictions obtained for low concentration values in some variables of interest. For example, negative predicted values were obtained for measured C content below 6.5gCkg⁻¹ soil, due to the wide range of concentrations used in the analysis (2.6–110.8gCkg⁻¹ soil). Several authors have been able to overcome this limitation, improving their predictions by using machine learning based non-linear regression methods, covariate adjustment (i.e., sand content) and by dividing the global library into small homogeneous datasets (Nocita et al., 2014; Rial et al., 2017; Ward et al., 2019).