A unified representation of the temperature dependences of soil microbial growth and respiration

Current temperature-dependence models

The Arrhenius model⁴⁰ is an empirical framework based on the van’t Hoff equation of equilibrium for chemical reactions. It was originally developed to describe how reaction rates increase with temperature, introducing the concept of activation energy. The model has since been widely applied across various scientific disciplines, including soil microbial ecology, to model temperature-dependent biological processes. The model assumes constant energy of activation and describes microbial rates (r) as:

The model has 2 calibration parameters: \(a\) and \(\mu\), which are the empirical Arrhenius factor and the activation energy or temperature characteristic of the reaction, respectively. In addition, \(R\) is the universal gas constant and \(T\) is the absolute temperature.

The Macro-Molecular Rate Theory (MMRT)⁴⁵ extends the Arrhenius’ theory further by incorporating the temperature dependence of heat capacity in large macromolecules, providing a more nuanced description of temperature effects on reaction rates. Originally developed to describe enzyme kinetics in pure cultures and biochemical systems⁴⁵, MMRT has been applied to soil microbial processes^16,46 and defines rates as:

The model has 4 calibration parameters, where \(\Delta {H}_{0}^{{{\ddagger}} }\), \(\Delta {C}_{{{{\rm{p}}}}}^{{{\ddagger}} }\) and \(\Delta {S}_{0}^{{{\ddagger}} }\) are the difference in enthalpy, heat capacity and entropy between the ground state and the transition state, respectively, and \({T}_{0}\) is the reference temperature (see refs. ^45,46 for details). \({k}_{{{{\rm{B}}}}}\) and \(h\) are the Boltzmann’s and Planck’s constants, respectively.

The MMRT-2S⁴⁹ model builds upon the original MMRT formulation (Eq. (2)) by introducing a temperature-dependent \(\Delta {C}_{{{{\rm{p}}}}}^{{{\ddagger}} }\):

where \(\Delta {C}_{{{{\rm{p}}}},{{{\rm{l}}}}}^{{{\ddagger}} }\) and \(\Delta {C}_{{{{\rm{p}}}},{{{\rm{h}}}}}^{{{\ddagger}} }\) are the heat capacities at the low- and high-temperature ranges, respectively, \(\Delta \Delta {H}^{{{\ddagger}} }\) is the enthalpy difference between the two ranges, and \({T}_{{{{\rm{c}}}}}\) is the temperature at the midpoint of the transition. MMRT-2S has 7 calibration parameters.

The Ratkowsky model¹² was initially defined as a linear relationship between the square root of rates and temperature, based on empirical observations⁴³:

This model was developed using microbial culture data and has been widely adopted to characterize microbial growth responses, particularly under suboptimal temperatures. It was later extended to describe the rates throughout the entire temperature range by introducing an exponential limiting factor:

This extended model involves 4 calibration parameters, where \(b\) and \(c\) are empirical coefficients, and \({T}_{\min }\) and \({T}_{\max }\) are the minimum and maximum temperatures at which microbial rates are conceptually assumed to be zero, respectively. Ratkowsky models have been increasingly applied^11,29,41,56 in soil microbial modeling due to their simplicity and robustness across temperature gradients. For improved estimates of these thermal parameters, authors recommended to first use Eq. (4) on rates at temperatures under optimal (e.g., <25 °C) to calibrate \(b\) and \({T}_{\min }\), and then use Eq. (5) for \(c\) and \({T}_{\max }\) on the whole range. Note that when temperature is much lower than \({T}_{\max }\), the contribution of the second half of the equation is negligible, and Eqs. (4) and (5) become identical.

The dual-kinetics Ratkowsky model

The new version of the Ratkowsky model defines rates as:

The estimate of the tipping point of the curves, representing the optimal temperature for growth and the temperature of slope acceleration for respiration, can be obtained by solving the derivative of Eq. (6a) with respect to temperature and setting it to zero:

It is noteworthy that the only difference between the proposed formulation and the original Ratkowsky model in Eq. (5) is in the sign preceding the exponential function in Eq. (6b). The Dual-Kinetics Ratkowsky model thus retains the 4 calibration parameters in the original Ratkowsky model (\(b\), \(c\), \({T}_{\min }\) and \({T}_{\max }\)). The negative sign (−) before the exponential function serves to account for the decreasing trend overserved in growth relationships beyond optimal temperatures. Conversely, the positive sign (+) for respiration is intended to capture the more pronounced slope in rates. The model thus integrates a biologically-plausible representation of anabolic and catabolic decoupling by capturing the adverse impact of high temperatures on microbial metabolism, functioning and viability, and how respiration is fueled by the increased availability of resources from destruction of living cells. Moreover, the model facilitates the estimation of the key thermal parameters \({T}_{\min }\), \({T}_{{{{\rm{tp}}}}}\), and \({T}_{\max }\), which are intended to characterize microbial trait adaptation to climate.

As in the original Ratkowsky model, two calibration approaches can be applied: a one-step calibration, which aims to achieve the best fit across the entire temperature range, and a two-step calibration (first rates <25 °C and then the whole range), which prioritizes the estimation of thermal parameters.

Experimental dataset for temperature relationships

To assess the fitness and biological relevance of the new formulation, we utilized a dataset sourced from Cruz-Paredes et al.¹¹, which included microbial growth (incorporating bacterial and fungal growth) and respiration temperature relationships from 70 soils. The surveyed soils encompassed a wide range of mean annual soil temperatures (from 2 °C to 20 °C), land uses (including forests, shrublands, grasslands, and croplands), soil organic matter contents (ranging from 3% to 94%), and pH levels (from 3.8 to 8.0), providing a broad representation of diverse microbial communities and soils. Rates were estimated under laboratory-controlled conditions at 5 °C intervals, spanning from 0 to 45 °C. Because these datasets were derived from short-term incubations of real soils, the temperature relationships used inherently reflect site-specific differences in microbial physiology and adaptation to different temperature regimes, but also factors in differences in resource availability and other intrinsic soil characteristics.

Temperature dependence models evaluation and comparison

Model parameters of Arrhenius, Ratkowsky, MMRT, MMRT-2S, and Ratkowsky DK were optimized by minimizing the sum of squared differences between observed and predicted values for each temperature dependence curve. To assess model performance, we estimated three key performance metrics—Coefficient of Determination (R²), Root-Mean-Square Error (RMSE), and Akaike Information Criterion (AIC)—for microbial growth, respiration, and the combined dataset (considering both microbial processes) (Table 1). Metrics were first computed separately, with rates normalized by their respective maximum to account for differences in magnitude across sites. For the joint evaluation of metrics for growth and respiration, normalized rates were combined, and the number of model parameters was summed to compute a combined AIC. These performance metrics were then averaged over the 70 sites, providing insights into the models’ goodness across a diverse range of soils. The Ratkowsky-based models were additionally evaluated using the two-step calibration approach to assess how prioritizing the estimation of thermal traits influenced model performance.

The Coefficient of Determination indicates the proportion of the variance in the observations that can be predicted, as:

where \({O}_{i}\) and \({P}_{i}\) are the observed and predicted rates at different temperatures, respectively, and \(\bar{{O}_{i}}\) is the mean of the rates in each site. A value closer to 1 indicates a better fit.

The normalized Root-Mean-Square Error quantifies model accuracy as the average magnitude of the difference between observations and predictions:

where \(N\) is the number of datapoints (rates at different temperatures). Lower RMSE values indicate better accuracy.

The Akaike Information Criterion balances goodness of fit with the number of parameters (\(p\)) calibrated for each model. The formula

is adjusted for the relatively small number of observations compared to parameters. Lower AIC values indicate a more parsimonious and effective model.

Modeling the impact on soil organic carbon stocks

To assess the potential impact of model selection on soil organic carbon stocks, we extended our comparative analysis between model outputs and observational datasets to now estimate how the cumulative fluxes of microbial respiration and growth influence soil carbon cycling in the field. To do this, we first modeled the annual cumulative fluxes at the 70 field study sites (in situ) by combining the intrinsic temperature relationships and soil temperature data at hourly time resolution. We then compared the results obtained by using the dependences derived from the models (Arrhenius, Ratkowsky, MMRT, MMRT-2S, and Ratkowsky DK) to those directly from the experimental temperature relationships (i.e., used smoothed data interpolation instead of one of the four models). We then estimated how the mismatch between models and raw data affected the equilibrium between microbial respiration (carbon released to the atmosphere as CO₂) and growth (microbial biomass retained as organic carbon within the soil). Finally, we employed a linear regression model to estimate the implications of model choice for soil carbon stocks along a gradient of mean annual soil temperature, based on the 70 sites, and mapped the potential effects across Europe at a 30-arc-second resolution based on current temperature and on a for future scenario of climate change predicted by the end of the twenty-first century.

To do this, first, hourly soil temperature (\({{{\rm{HST}}}}\)) was modeled at each site using the equation:

where soil temperature oscillates around the mean annual soil temperature (\({{{\rm{MAST}}}}\)), with the second and third terms of the equation representing seasonal and diurnal variation, respectively. \({A}_{{{{\rm{s}}}}}\) and \({A}_{{{{\rm{d}}}}}\) describe the amplitude of this seasonal and diurnal fluctuations, while \({t}_{{{{\rm{d}}}}}\) is the time in days and \({t}_{{{{\rm{h}}}}}\) is the time in hours. The values for \({{{\rm{MAST}}}}\), \({A}_{{{{\rm{s}}}}}\) and \({A}_{{{{\rm{d}}}}}\) for each site were obtained from the “SoilTemp” database, a global 30-arc-second resolution soil temperature dataset by ref. ⁶⁹. that provides estimates in the top 0–5 cm of soil. The seasonal (\({{{{\rm{lag}}}}}_{{{{\rm{s}}}}}\)) and diurnal (\({{{{\rm{lag}}}}}_{{{{\rm{d}}}}}\)) lag phases were set to 243 days and 14 h, respectively, to simulate the natural cycles in soil temperature.

We then employed a mass-balance conceptual approach that quantified how the accuracy of the temperature-dependent models potentially affects the key fluxes driving changes in soil carbon stocks. This approach assumes that soil carbon stocks at each location are in a state of pseudo-equilibrium, meaning that annual changes are minimal, with the total carbon released from the soil through microbial respiration balanced by the carbon incorporated into the soil through microbial growth. The equilibrium was expressed as the ratio:

where \(\alpha\) represent the quotient of equilibrium, reflecting the proportionality between the annual cumulative fluxes for respiration (\({\left.{{{\rm{cumR}}}}\right|}^{{{{\rm{d}}}}}\)) and growth (\({\left.{{{\rm{cumMG}}}}\right|}^{{{{\rm{d}}}}}\)) obtained using a locally weighted regression technique on temperature-relationship empirical data to ensure smooth transitions and reduce noise. While this interpolation did not fully capture the system’s nonlinearity, it allowed for a standardized comparison between models and raw data, similarly to how model fitness was analyzed.

We then estimated how each model deviated from the equilibrium established using the raw data with the following expression:

where \({\left.{{{\rm{cumR}}}}\right|}^{{{{\rm{m}}}}}\) and \({\left.{{{\rm{cumMG}}}}\right|}^{{{{\rm{m}}}}}\) are the annual cumulative fluxes for respiration and growth obtained from the models. It is important to note that these estimates of the potential effects on soil organic carbon stocks should not be interpreted in absolute terms, as they represent a comparison of how the accuracy of each model may influence the stocks. In addition, the estimates focused solely on temperature effects, as other environmental factors such as moisture content or oxygen availability were not considered. Our approach thus compares model predictions to a reference based on temperature relationships derived directly from experimental data. Although these relationships are based on precise laboratory measurements, small experimental errors are amplified when propagated into SOC stock estimates. As a result, the deviation values (\({{{\rm{Dev}}}}\)) may overstate the true uncertainty and should be interpreted as a relative measure of model performance rather than as definitive estimates of absolute error. Values closer to zero indicate a better alignment of the model with the raw data, predicting fluxes that are closer to the equilibrium between carbon emissions and carbon sequestration. Positive values for \({{{\rm{Dev}}}}\) suggest that the model is overestimating carbon emissions, while negative values indicate that the model is underestimating them (or overestimating carbon buildup in soil).

Finally, we employed a linear regression analysis to examine how mean annual soil temperature influenced the deviation calculated from the five different temperature-dependence models. By analyzing the relationship between them across the study sites, we generated a map of the projected effects of model choice on potential fluxes altering soil organic carbon stocks across Europe, based exclusively on the mean annual soil temperature of each site and without applying regional smoothing. This map was produced at a high spatial resolution of 30 arc seconds, based on current soil temperature data (‘current climate’) and projected future climate change scenarios (‘future climate’). For the future climate scenario, we considered an estimated increase in global average soil temperature of +5.8 °C at a depth of 0–5 cm (3.5 cm) by the end of the twenty-first century. This estimate was derived from an EC-Earth3-Veg model simulation under the high-emissions scenario SSP585, with soil temperature data retrieved from the Earth System Grid Federation (ESGF) repository (https://aims2.llnl.gov/search/cmip6/, accessed 19 August 2024). The projection enabled us to assess how changes in temperature regimes might affect the model-driven discrepancies in soil carbon flux estimates, offering insights into the potential effects of temperature-dependence model choice under different climate conditions, for instance, in areas warmer than Europe. This estimation did not account for further adaptation of microbial traits to higher temperatures but instead used current thermal traits observed along the European gradient.

Statistics

The Kruskal–Wallis test was employed to assess differences in the key metrics for model performance (R², RMSE and AIC), enabling the identification of the models that provided the best results in terms of goodness of fit and parsimony. Subsequently, a pairwise post hoc analysis with Dunn-Sidak adjustment was conducted to compare models and group them based on statistically significant differences in performance (p > 0.05). The best methods within each category were identified and classified as group ‘a’ (superscript letter). To test the correlation between the parameters obtained from each model and the mean annual soil temperature, a Spearman correlation analysis was performed. In addition to the p value and correlation coefficient (\(\rho\)), we calculated a pattern score (\(-\log \left(p\right)\cdot {\rho }^{2}\)) to capture both the strength and reliability of the trend. Subsequently, linear regression was used to evaluate the relationship between \({T}_{\min }\), \({T}_{{{{\rm{tp}}}}}\), and \({T}_{\max }\) of the Ratkowsky DK model and temperature across the gradient.

All calculations and statistical analyses were performed using MATLAB (R2024b).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.