General comments:

Based on national scale soil samplings, this manuscript proved the potential of implementing mid-IR spectra and machine-learning for MAOC and C deficit prediction. The results show that the CUBIST models for both MAOC and C deficit prediction have good performance, advocating their future application. They also make these models interpretable by matching absorption features of the mid-IR spectra and coefficients in models among different modeling rules. Nevertheless, several issues raised during my review which I think should be addressed before publication.

• The investigation of model interpretability should be modified. Since the SHAP values coincide with the regression coefficients of the CUBIST rules, there is large redundancy between the SHAP analysis and that of CUBIST rules demonstration. In other words, the interpretation that positive SHAP values had a positive impact on the model prediction also applies to that of coefficient values in multivariate regression. The authors should demonstrate the additive value of the SHAP analysis. In addition, if the authors manage to do so, then they should also perform the SHAP analysis on MAOC prediction model. Otherwise, the authors should declare the reason why they only perform the SHAP analysis on C deficit prediction model. In addition, the so-called interpretability stops by pointing out impactful wavenumber and its chemical identity. The interpretability should involve more explanatory descriptions. For instances, in line 259, “absorptions for quartz and other minerals in the fingerprint region were also important in the models, but negatively affected the estimates”. What did this result tell us? Is that because the relatively larger amount of quartz likely indicates a sandy texture of soils, thus indicating less mineral capacity and likely low C deficit?

• The discussion section should be modified in several aspects. First, the authors stated that the spectroscopic approach enables many more measurements than conventional methods, enhancing our understanding of how MAOC and C deficit vary in the soil in space and time. However, the approach that this study implemented still involved destructive samplings over large geographical scale, which still belong to conventional methods. In other words, in order to monitor C deficit dynamics over time, researchers need long-term large-scale samplings to get the new mid-IR spectra from soils, even they have built the CUBIST models. Therefore, the statement will be a better fit for spectroscopic approaches which use spectra from non-destructive remote sensing techniques, i.e. spectra from satellites, even though the model accuracy of these studies tends to be lower than this study. If insist using the statement mentioned above, the authors should point out the potential that laboratory-based spectroscopic approaches can help improve the performance of that of remote sensing spectroscopic approaches. Second, the authors pointed out that the frontier line approach can have a more accurate estimate of MAOC maximum capacity than that of quantile regression in discussion part. However, Shi et al (doi.org/10.1016/j.geoderma.2025.117181) has implemented a local approach for the quantile regression method, which has the merit of avoiding under- or over-estimations. The authors should incorporate Shi’s study into the discussion section and modify the relevant statements.

• The particle size of clay and silt content and of fine fraction in soil fractionation are methodologically mismatched, which induced errors. The mineral capacity between soil particles under 20 μm and 50 μm are different. Because these two sets of soil minerals have different structure in their components. For instance, the 50 μm set might constitute more quartz, feldspar, and 1:1 type clay mineral, which have lower C absorption capacity than that of 2:1 type clay mineral. Thus, the C absorption capacity of soil minerals partitioned by 20 or 50 micrometers cannot represent each other. Using 20 μm clay and silt content to capture MAOC maximum capacity corresponding to 50 μm fractionation protocols does not robustly reflect the relationship. There might be a few options for improvement: changing the model for clay and silt prediction, laboratory work for clay and silt content, or at least acknowledging this limitation in the discussion.

Minor comments:

Line 41: Instead of fitting 90^th quantile regression, Georgiou et al used 95^th quantile regression. Please check.

Line 116: Did this back-transformation be performed during uncertainty analysis? Since the authors used logarithm when fitting the frontier line, the upper and lower uncertainty intervals would be different between that undergone first calculating intervals then back-transformation, and that undergone first back-transformation then calculating intervals. Please clarify.

Line 124: What specific are the offset corrections? SNV transformation is well-known in spectroscopic area, while offset correction tend to be a series of mathematical operation on the spectra. Please clarify or at least provide reference.

Line 174-176: The result is not intuitive. It is hard to tell whether samples in Rule 3 have higher absorption in the 2946–2850 cm⁻¹ region than that of Rule 4, given the scale of the y-axis in the two plots are not consistent. Could the authors please make this comparison more intuitive, thus better supporting the statement?

Line 255: The authors mentioned they have propagated the uncertainties from the frontier lines fits and the CUBIST models to our final predictions. Do the uncertainties of the frontier line fits have anything to do with the uncertainty of C deficit CUBIST model? Because the latter is demonstrated with parameters like RMSE only for C deficit model not its upper or lower 95% confidence intervals CUBIST models. There is a mismatch between the grey areas in Figure 5 and statistical parameters of the C deficit CUBIST model, indicating there is no propagation of the intervals to the final C deficit prediction. Please clarify.

The manuscript is clearly written. Analyses are properly conducted. I have just a few questions.

1. Since clay and silt content are also predicted from MIR spectra, how does the accuracy of these predictions influence the calculation of Cdef and spectral modeling?
2. The caption of Figure 3 – the last sentence is repeated.
3. What is the direct linear or nonlinear relationship between MAOC and Cdef?
4. Since many spectral regions are identified that relate to organic groups, clay, and quartz, what is the model accuracy when using these soil properties to directly predict Cdef?
5. Since leave-site-out cross-validation is used in the study, how does the model accuracy compare when an independent validation is applied?
6. How does model performance vary with the three depth layers?