**Baseline microbe-enzyme model representing explicit resource allocation trade-off**

Our baseline microbe-enzyme model (fig. S10) is a second-order model of soil organic carbon (SOC) dynamics that expands upon the model of Allison et al 2010 (AWB) by explicitly representing a microbial trade-off between carbon allocation to growth versus extracellular enzyme production. It includes SOC (C), dissolved organic C (DOC), microbial biomass C (M) and explicit enzyme mass C (Z). SOC increases with litter input (*I*), and declines with leaching (*e**C*) and enzyme decomposition, which is assumed to be catalyzed according to Michaelis-Menten kinetics of the enzyme pool as represented by the last term in Eq. 1:

(1)\(\frac{dC}{dt}=I-{e}_{C}C-\frac{{v}_{max}^{D}C}{{K}_{m}^{D}+C}Z\)

DOC receives the product of SOC decomposition, as well as the product from the recycling of dead microbes (*d**M*) and inactive enzymes (*d**Z*). We did not include litter input to the DOC pool, represented in AWB, because the steady states are not sensitive to it (see sensitivity analysis in Table S2) and it simplifies the analytic solution. DOC decreases with leaching (*e**D*) and microbial uptake, where it is assumed that microbial cells have a limited amount of transporters, hence another Michaelis-Menten process represented by the last term in Eq. 2:

(2)\(\frac{dD}{dt}=\frac{{v}_{max}^{D}C}{{K}_{m}^{D}+C}Z+{d}_{M}M+{d}_{Z}Z-{e}_{D}D-\frac{{v}_{max}^{U}D}{{K}_{m}^{U}+D}M\)

Microbial biomass turns over with rate *d**M*, and grows proportionally to DOC uptake and a growth efficiency term that accounts for the loss to respiration due to the energetic cost of biomass production (*γ**M*), minus a portion of it allocated to enzyme production (*φ*):

(3)\(\frac{dM}{dt}=\left(1-\phi \right){\gamma }_{M}\frac{{v}_{max}^{U}D}{{K}_{m}^{U}+D}M-{d}_{M}M\)

Finally, the enzyme pool grows with microbial uptake and allocation to enzyme production, and declines with enzyme turnover rate:

(4)\(\frac{dZ}{dt}=\phi {\gamma }_{Z}\frac{{v}_{max}^{U}D}{{K}_{m}^{U}+D}M-{d}_{Z}Z\)

The resource allocation trade-off is parameterized with *φ*, with *φ* the fraction of DOC uptake allocated to enzyme production, and (1-*φ*) the fraction allocated to yield. Parameter units and values are detailed in Table S1.

**Steady states**

The model (Eqs. 1–4) possesses either one globally stable equilibrium, or three equilibria (one of which is always unstable). There are thresholds *φ*min and *φ*max such that the single globally stable equilibrium exists for *φ* < *φ*min or *φ* > *φ*max and is given by *C* = *I*/*e**C*, *D* = 0, *M* = 0, *Z* = 0. Thus, at this equilibrium, the microbial population is absent and no decomposition occurs. For *φ*min < *φ* < *φ*max, the microbial population can either go extinct (then the system stabilizes at the same equilibrium as before) or persists at or around a non-trivial equilibrium, which can be solved analytically (more details on the non-trivial equilibrium, *φ*min and *φ*max in Note S1).

**Temperature dependence**

At the molecular and cellular level, the effect of warming on microbial decomposition is mediated by the temperature sensitivity of intra- and extra-cellular enzymatic activity (*7*, *18*, *32*). In the baseline ‘kinetics-only’ scenario, microbial uptake parameters (maximum uptake rate, *v**U**max*, half-saturation constant, *K**U**m*) and exoenzyme kinetics parameters (maximum decomposition rate, *v**D**max*, half-saturation constant, *K**D**m*) increase with temperature (*33*, *34*) in an exponential manner:

(5)\({v}_{max}^{D}= {v}_{0}^{D}{e}^{-\frac{{E}_{v}^{D}}{R\left(T+273\right)}}\)

(6)\({K}_{m}^{D}= {K}_{0}^{D}{e}^{-\frac{{E}_{K}^{D}}{R\left(T+273\right)}}\)

(7)\({v}_{max}^{U}= {v}_{0}^{U}{e}^{-\frac{{E}_{v}^{U}}{R\left(T+273\right)}}\)

(8)\({K}_{m}^{U}= {K}_{0}^{U}{e}^{-\frac{{E}_{K}^{U}}{R\left(T+273\right)}}\)

**Eco-evolutionary optimization of microbial resource allocation**

Assuming heritable variation in the exoenzyme allocation fraction trait, *φ*, we used the framework of adaptive dynamics (*35*, *36*) to predict the strength and direction of selection on trait *φ* and the evolutionarily stable value, *φ**. In this framework, eco-evolution is modeled as a competition process between a ‘resident strategy’ (wild-type) and alternate strategies (‘mutants’) within a set of feasible phenotypes. In a given environment (*e.g.* at a given temperature), an evolutionarily stable strategy (ESS) is a phenotype that, when resident, no mutant can invade. The adaptive dynamics framework provides the mathematical criteria to identify ESSs and check their attractivity, i.e. that they can be reached by a sequence of small evolutionary steps, each step involving the replacement of a resident phenotype by a mutant phenotype. Here the set of feasible phenotypes is the range (*φ*min, *φ*max) at a given temperature, for which the non-trivial ecosystem equilibrium exists.

To model the competition effect of a resident phenotype, *φ*res, on the population growth of a mutant phenotype, *φ*mut, we extended the baseline microbial-enzyme model written for a single type (Eq. 3). To account for the local nature of the interaction between rare mutant and common resident cells, we introduced a function (hereafter denoted by *c*) of the difference between *φ*res and *φ*mut to measure how local decomposition by mutant and resident cells differ from ‘mean field’ (average) decomposition by resident cells. Thus, for given *C*, *D*, *Z*, the growth of the mutant population is governed by:

(9)\(\frac{d{M}_{mut}}{dt}=\left(1-\phi \right){\gamma }_{M}\frac{{v}_{max}^{U}\left(1+c\left({\phi }_{mut}-{\phi }_{res}\right)\right){D}_{res}}{{K}_{m}^{U}+\left(1+c\left({\phi }_{mut}-{\phi }_{res}\right)\right){D}_{res}}{M}_{mut}-{d}_{M}{M}_{mut}\)

where *D*res is the equilibrium *D* predicted by the ecosystem model for the resident phenotype *φ*res. Here function *c* satisfies *c*(0) = 0, *c*(*z*) > 0 if *z* > 0 and *c*(*z*) < 0 if *z* < 0. The underlying assumption is that each microbe has access to DOC partly as a public good and partly as a private good (*37*). The public good part results from the diffusion of exoenzymes. The private good part results from local decomposition at the microscopic scale of cells and exoenzymes that they produce themselves. A mutant cell that invests more in exoenzymes has access to more DOC than the average resident cell because the cell’s private good is greater whereas all cells share the same public good. In a spatially implicit model like ours, diffusion is not directly modeled, but its effect on the accessibility of DOC to a mutant strain can be phenomenologically accounted for by a parameterization that puts mutant cells at a competitive advantage for DOC if the mutant phenotype invests more in exoenzyme production than the resident phenotype, or at a competitive disadvantage if the mutant phenotype invests less. This parameterization is achieved with the function *c* in Eq. 9, where *c* < 1 when *φ*mut < *φ*res and *c* > 1 when *φ*mut > *φ*res. This phenomenological approach is consistent with the mathematical construction and numerical analysis of a spatially explicit model of resident-mutant local interaction that accounts for soil diffusion (*38*).

Mutant relative fitness *s*(*φ**mut* - *φ**res*) is given by the mutant population growth rate per unit biomass:

(10)\(s\left({\phi }_{mut},{\phi }_{res}\right)=\left(1-{\phi }_{mut}\right){\gamma }_{M}\frac{{v}_{max}^{U}\left(1+c\left({\phi }_{mut}-{\phi }_{res}\right)\right){D}_{res}}{{K}_{m}^{U}+\left(1+c\left({\phi }_{mut}-{\phi }_{res}\right)\right){D}_{res}}-{d}_{M}\)

The selection gradient is then obtained by taking the first order derivative of the invasion fitness with respect to the mutant trait:

(11)\(\nabla s\left(\phi \right)=\frac{{d}_{M}}{1-\phi }\left(\left(1-\phi -\frac{{d}_{M}}{{v}_{max}^{U}{\gamma }_{M}}\right){c}_{0}-1\right)\)

where *c*0 = *c’*(0) measures the local competitive advantage to stronger exoenzyme producers, that we call ‘competition asymmetry’. Note that by definition of function *c*, we always have *c*0 > 0. Variation in *c*0 may be caused by different soil diffusion properties, due to *e.g.* physical texture or moisture.

Trait values that nullify the selection gradient are called ‘evolutionary singularities’. An evolutionary singularity can be attractive or repelling, and invadable or non-invadable. Evolutionary singularities that are attractive and non-invadable represent potential end-points of eco-evolution. Evolutionary singularities that are attractive and invadable can lead to evolutionary branching (*36*). In a given environment (fixed parameters, constant temperature) there is at most one evolutionary singularity given by defining *φ** as the value of *φ* that makes \(\nabla s\left(\phi \right)\)= 0:

(12)\({\phi }^{*}=1-\frac{{d}_{M}}{{\gamma }_{M }{v}_{max}^{U}}-\frac{1}{{c}_{0}}\)

Existence of *φ** > 0 requires \(\frac{{d}_{M}}{{v}_{max}^{U}{\gamma }_{M}}<1\) and \({c}_{0}>\frac{1}{\left(1-\frac{{d}_{M}}{{v}_{max}^{U}{\gamma }_{M}}\right)}\). Thus, the (cooperative) trait *φ* can evolve above zero only if the local competition advantage to stronger enzyme producers is large enough. The condition for *φ** to be evolutionarily stable is *c’’*(0) < 2 c02, which means biologically that the competition advantage needs to be towards lowering the gain for cheaters rather than increasing it for producers. No other condition than existence is required for *φ** to be always convergent. Here we assumed that function *c* is such that *φ** is evolutionarily stable and attractive.

**Calculations at the global scale**

We obtained global projections of soil carbon stocks by calculating the steady states on a grid of 1° × 1° across terrestrial continents. Projected soil temperature was obtained from the version 4 of the Community Climate System Model (CCSM4) for the Representative Concentration Pathway 8.5 (RCP8.5) from 2005 to 2100 (publically available online at https://www.cesm.ucar.edu/models/ccsm4.0/) (*11*). It provides daily predictions across 192 points of latitude and 288 points of longitude. We averaged those projections over one year each decade to calculate the local steady states used in all figures. For all sites where the microbial activity is positive, the soil C stock is estimated as the sum of the steady states of SOC and microbial biomass (in units of carbon mass per m2 per centimeter of soil depth). To calculate the effect of optimization, we calculated the sum of changes in soil C stock of every decade between 2010 and 2100. For the predictions obtained with alternative scenarios of microbial temperature sensitivity, we used equations and parameters from Allison et al. (*5*) and Hagerty et al. (*6*) (see equations in Note S2). For the predictions obtained with spatially and temporally variability in litter input, we used the litter projections from the CESM (same model source as for temperature), and averaged them as well over a year each decade in each site. Finally for the predictions obtained with spatially variable temperature sensitivity of enzyme kinetics, we used the 5 sets of sensitivity parameters obtained by German et al. (*7*). We divided continents into 5 clusters based on temperature and moisture (fig. S8), to which one of the 5 sets of enzyme temperature sensitivity parameters was attributed (Table S3).