Metabolic rates of prokaryotic microbes may inevitably rise with global warming

Understanding how the metabolic rates of prokaryotes respond to temperature is fun-damental to our understanding of how ecosystem functioning will be altered by climate change, as these micro-organisms are major contributors to global carbon efflux. Ecological metabolic theory suggests that species living at higher temperatures evolve higher growth rates than those in cooler niches due to thermodynamic constraints. Here, using a global prokaryotic dataset, we find that maximal growth rate at thermal optimum increases with temperature for mesophiles (temperature optima ≲ 45°C), but not thermophiles (≳ 45°C). Furthermore, short-term (within-day) thermal responses of prokaryotic metabolic rates are typically more sensitive to warming than those of eukaryotes. Given that climatic warming will mostly impact ecosystems in the mesophilic temperature range, we conclude that as microbial communities adapt to higher temperatures, their metabolic rates and therefore, carbon efflux, will inevitably rise. Using a mathematical model, we illustrate the potential global impacts of these findings.

at different time-scales on ecosystem fluxes studied. Figure 1: Three alternative hypotheses for the short-vs. long-term responses of thermal performance curves of a fitness-related metabolic trait in response to environmental warming. A. Hotter-is-Better: organisms adapt around a global, inter-specific, thermal constraint (black line, Boltzmann-Arrhenius fitted to intra-specific curve peaks), such that the average intra-specific (short-term) activation energy (Ē S ) is statistically indistinguishable from the inter-specific (long-term) activation energy of the group of organisms (E L ), and both are greater than zero. See methods for more details on the the definition and estimation ofĒ S and E L , and the statistical methods used to differentiate between them. Note that each intra-specific TPC represents the short-term thermal response of each organism's population. Inset panel illustrates how this would look in an Arrhenius plot. B. Equalisation of fitness: selection overrides thermodynamic constraints, such that trait performance at T pk is on average the same (E L = 0). Alternatively the same effect of E L = 0 may occur due to or thermodynamic constraints on enzymes in fact restricting metabolic rate (and therefore fitness) at higher temperatures. C. Weak biochemical adaptation: an intermediate scenario where E L > 0 but significantly less thanĒ S . Panel C also illustrates the the Sharpe-Schoolfield TPC model parameters (eqn. 1, Methods).
Here we build and analyse a global dataset of TPCs in bacteria and archaea to quantify general patterns 69 in both the short-term (intra-specific) response, and to test whether the HiB hypothesis holds (long-term, 70 inter-specific response) within and across taxonomic and functional groups adapted to different tempera- Fig. 1) at the same rate as r max would increase with temperature (parameter E S ), on average, within 92 single strain's TPC. Our analysis relies on P pk -T pk pairs across strains because data within strains are 93 largely lacking, and the HiB pattern is expected to apply across strains within monophyletic taxonomic 94 groups (such as archaea and bacteria) 24,30 . Analysing this relationship across 416 bacterial and 82 95 archaeal strains, we find that hotter is indeed better (HiB holds) across mesophiles (Ē S and E L are > 0, 96 and their 95% CIs overlap; Fig 2 and Table 1). However, this result does not extend to thermophiles, 97 where instead fitness is on average invariant with respect to temperature. Thermophiles have evolved 98 specific adaptations to extreme temperature stress, such as mechanisms to cope with increased membrane 99 permeability at high temperatures 31 and thus adaptation to such niches may incur a fitness cost to 100 thermophiles as seen in our results. This result is in concurrence with an investigation of the maximum 101 growth rates of life on Earth, which found increases in microbial growth up to a peak before an attenuation 102 of growth rates in warmer adapted organisms 28,29 . Table 1: Estimated mean short-term and long-term activation energies across archaea and bacteria, and test of the HiB hypothesis. Estimated mean E S and E L values (95% CI ranges in parentheses) for bacteria and archaea split by thermal niche (also see Fig 2). The last column indicates whether or not the HiB hypothesis is supported.

Kingdom Thermal
Niche  Figure 2: Patterns of short-and long-term thermal responses of growth rate (fitness) for archaea and bacteria. A and B: Activation energies (with 95% Confidence Intervals) from Boltzmann-Arrhenius model fits (E L , blue triangle) compared to mean activation energy (Ē S , red square) of the intraspecific (short-term) thermal responses. Distribution of all the E S s is also shown (orange points). C and D: Arrhenius plots (x-axes inverted to aid visualisation) fitted to mesophile and thermophile sub-groups separately within bacteria and archaea, respectively. That is, the lines (the long-term thermal responses) are the Boltzmann-Arrhenius model fitted using weighted regression to mesophile and thermophile data separately, after determining the breakpoint (Methods). The HiB hypothesis is best-supported for the mesophile sub-group in both panels, while equalisation of fitness is best supported in the the thermophile sub-groups (Table 1).

105
Under the MTE, the global (inter-specific) thermodynamic constraint is expected to center around 106 0.65eV 5,6 . Mean intra-specific thermal sensitivities have been found to be very similar to this value, 107 although the distribution is right-skewed with a median value of ∼0.55eV 2 . In contrast, we find that 108 mean thermal sensitivities for both bacteria and archaea fall significantly above 0.65eV (Fig. 4, bac-109 teriaĒ S = 0.88eV; archaeaĒ S = 0.95eV). We observe the same right-skew in activation energies for 110 prokaryotes as seen across other organisms and traits 2 . Even after accounting for this skew by taking the 111 median instead of the mean, activation energy still falls significantly above 0.65eV (bacteria median = 112 0.84eV, archaea median = 0.80eV; Fig. S2). Furthermore, we see a consistent pattern of median thermal 113 sensitivity >0.65eV throughout the lower taxonomic groupings (Fig. 3A).

114
To further understand these findings we also categorised the data into various groups based on func-115 tional traits (Fig. 3B). Again we find mean and median thermal sensitivity >0.65eV in the majority of 116 functional groups, suggesting that this high E is a trait generally conserved across prokaryotic organisms.
117 Figure 3: Variation in thermal sensitivity among prokaryotic groups. Comparison of intraspecific population growth rate activation energies (Ē S ) across taxonomic levels (A) and functional trait groupings (B). Points and error bars represent weighted mean and 95% CIs of E S for each group. Groups shown are those with at least five data points, the number in brackets indicates the number of data points from whichĒ S was calculated for each grouping. The dotted line marks 0.65eV, the mean E previously reported within the MTE framework. Grey triangles mark the median E S for each group.
Here we have focused on the TPCs (and activation energies) of population growth rate. However, to Potential Ecosystem-level Impacts

130
Our results higher sensitivity of both short-(higher intraspecific activation energies) and long-term (higher  depending on the quantity of prokaryotic biomass already in the system. 156 Figure 5: Potential changes in climate-driven short-and long-term ecosystem carbon flux due to difference in sensitivity between prokaryotic and eukaryotic thermal responses. A. Heat map of % short-term increase in flux with 10 • C temperature increase of model ecosystems with bacteria having a different activation energy on average than eukaryotes, relative to ecosystems with all components having the same (0.65eV) average activation energy. The flux change is shown over a range of ecosystem biomass compositions in terms of heterotrophs vs. autotrophs and bacterial proportion of the heterotrophs. The scale of emergent activation energies and Q 10 s for the ecosystems with amplified flux are also shown. B. Similar to A, but for long-term flux increase under a 4 • C warming scenario. Values for the short-and long-term thermal sensitivity of bacterial thermal responses used in these calculations are our estimatedĒ S and E L respectively for mesophilic bacteria ( Table 1). The mathematical model is described in Methods.

Discussion
is-Better constraint which results in the flux at thermal optimum also increasing with (longer-term) 164 adaptation. Moreover, this coupling across timescales is expected not just in methanogens, but across 165 most major mesophilic prokaryotes, including those involved in aerobic respiration. HiB curve (Fig. 1A). Alternatively, species sorting may occur such that prokaryotes inherently better-171 adapted to higher temperatures take advantage of temperature increases. This would have the same 172 overall effect because these prokaryotes would also effectively be further up the inter-specific temperature 173 response curve (Fig. 2). In either case, under HiB, we can expect global warming to result in prokaryotic 174 communities with higher metabolic rates on average. Thus overall, our results suggest that further 175 production of greenhouse gases from the prokaryotic component of ecosystems is likely to increase in 176 general, and at a greater rate than that by component eukaryotic organisms (Fig 5).

177
While in general, we see a tendency towards high thermal sensitivity (E S ) in prokaryotes, there are 178 taxonomic subgroups within our dataset for which this is not the case (Fig. 3). For example,Ē S for 179 mesophilic archaea as a whole does not deviate significantly from the MTE 0.65eV average (Table 1). This  To each TPC in the dataset, we fitted a modified Sharpe-Schoolfield model 52 (eq. 1): Here, T is temperature in Kelvin (K), B is a biological rate, B 0 is a temperature-independent metabolic 256 rate constant approximated at some (low) reference temperature T ref , E is the activation energy in electron 257 volts (eV) (a measure of "thermal sensitivity"), k is the Boltzmann constant (8.617 × 10 −5 eV K −1 ), T pk 258 is the the temperature where the rate peaks, and E D the deactivation energy, which determines the rate 259 of decline in the biological rate beyond T pk . We fit this model to individual TPCs and solve for T = T pk 260 to calculate the population growth rate at T pk (P pk ) for each strain. Note that this has been reformulated 261 from the model presented in the original paper, to include T pk as an explicit parameter 53 .

262
Each strain's TPC has a potentially different T pk and P pk . Compiling these values across strains 263 yields an inter-specific thermal response curve (Fig 1). We fit the Boltzmann-Arrhenius equation (eq. 2, 264 essentially the numerator in eq. 1) to these peak values to calculate inter-specific activation energy. NumPy package, using a least squares regression method to minimise the fits.

268
Comparing Short-and Long-term thermal responses 269 We determined whether Hotter-is-Better by testing whether the activation energies from intra-(short-270 term) and inter-specific (long-term) TPCs were (statistically) significantly different. For each intra-271 specific curve we fitted the Sharpe-Schoolfield model (eq. 1) and extracted the intra-specific activation 272 energy (E S ), peak temperature (T pk ) and corresponding growth rate (P pk ) for each curve. TPCs without 273 a peak are thus excluded from this analysis.

274
To estimate E L we fitted the Boltzmann-Arrhenius model (eq. 2) to the T pk and P pk values estimated 275 from the intra-specific thermal responses. To account for uncertainty in the original Sharpe-Schoolfield 276 model fits to the intra-specific curves, we fitted Boltzmann-Arrhenius using a weighted regression (see 277 accounting for uncertainty). In order to provide a comparison between intra-and inter-specific activation 278 responses, we used bootstrapping to generate confidence intervals (CIs) around the mean in each case.

279
To provide boostrapped CIs for E L from the modified Boltzmann-Arrhenius fits, the data was re-sampled  We then determined whether the data was consistent with either of the three hypotheses (main text 285 Fig. 1 S and E L don't overlap). Under a HiB scenario, P pk will increase with T pk across strains, and according other, but not zero. Alternatively, if growth rates are not constrained by thermodynamics and P pk does 294 not increase with temperature, then E L will be close to zero (CI for E L includes zero), and HiB can be 295 rejected. Finally, in scenarios where thermodynamic constraints may be partially evident but somewhat 296 overcome by adaptation,Ē S and E L will both be positive, but withĒ S being significantly greater than

299
Weighted means were used to account for uncertainty in Sharpe-Schoolfield point estimates when calculat-300 ingĒ S and when fitting inter-specific Boltzmann-Arrhenius curves. After performing Sharpe-Schoolfield 301 fits we extracted the E and µ pk point estimates as well as the covariance matrix. We then sampled 302 1,000 times from a bivariate distribution accounting for the covariance, producing 1,000 model parameter 303 combinations. We used these parameters to generate 1,000 different Sharpe-Schoolfield curves, providing 304 a distribution of E and µ pk from which we took the standard deviations (SD E and SD µ ) as a measure 305 of uncertainty. In some cases the Sharpe-Schoolfield fit did not produce a covariance matrix and these 306 fits were excluded from further analysis.

307
When combining E values across strains to calculateĒ S we, took weighted arithmetic means of E 308 to account for uncertainty in the original fits, where W eight = 1/(SD E + 1). Similarly, when fitting 309 the HiB hypothesis is accepted or not for different groupings, however we felt that it was important to 312 acknowledge and account for error in the underlying Schoolfield fits so that our results were not skewed 313 by poor parameter estimates from questionable fits, hence this step was included. Figure S2 illustrates 314 the differences betweenĒ S calculated with and without a weighting -applying a weighting pushesĒ S 315 down a little, likely due to high E values obtained from fits to lower quality data. In either case, with 316 or without a weighting,Ē S falls significantly above the 0.65eV MTE average activation energy for both 317 Bacteria and Archaea.

318
Taxonomic and Physiological Groupings

319
Pychrophiles and mesophiles inhabit low to medium temperature ranges, while thermophiles and hyper-320 thermophiles grow at much higher temperatures 54 . The distinction between mesophiles and thermophiles 321 is usually defined relatively arbitrarily, with mesophiles often considered strains with thermal optima up 322 to 45°C and thermophiles those with thermal optima of 55°C and above 54 . Corkrey et al. 28 found a peak 323 in microbial growth rates at ∼42°C (mesophile peak) followed by an attenuation of maximum growth rates 324 until a second peak at ∼67°C (thermophile peak), suggesting a biological transition between mesophiles 325 and thermophiles.

326
In order to determine whether it was appropriate to consider mesophiles and thermophiles separately, 327 we performed a break-point analysis on our dataset using the 'Segmented' R package 55 . Segmented is 328 not compatible with non-linear least-squares (nls) fitting, so this was performed with a linearised version 329 of Boltzmann-Arrhenius, i.e. x ∼ y where x = 1/(kT pk ) and y = log(µ pk ). As this process was merely 330 to confirm whether it was appropriate to split the data into mesophiles and thermophiles as suggested 331 by eye, it is not important that these linearised fits may give slightly different slope and intercepts to 332 the weighted nls fits. Using this methodology we determined significant break-points for bacteria and 333 archaea within our growth rates dataset at 40.48°C and 46.21°C respectively. These are similar to the 334 ∼42°C mesophile growth rate peak seen by Corkrey et al. 28 and were thus used as cut-off points for 335 defining mesophiles and thermophiles in our analysis.

336
In addition, archaea are typified by their adaptations to energetically demanding niches, while in 337 contrast bacteria perform better in more "ambient" environments 43 . A major physiological difference 338 between these taxa lies in their fundamentally divergent membrane structures. This affects these organ-339 isms' abilities to maintain proton gradients and thus drive metabolism under different conditions 43 , a 340 difference that may be particularly important for thermal performance. As such, we separate bacteria 341 and archaea in our analysis as disparate organisms with divergent evolutionary histories.

342
In order to classify prokaryotes by the energy generating metabolic processes that they use, we took note 343 of the growth conditions used when initially digitizing the TPC data. For the majority of heterotrophic 344 bacteria and archaea this was simply whether they were grown under aerobic or anaerobic (fermentative) 345 conditions. However there are also a number of strains utilising more exotic metabolic processes such as 346 methanogenesis, sulfur reduction, etc. In these cases we matched taxa against those able to utilise certain 347 metabolic reactions according to Amend & Shock 56 before manually checking the culture conditions in 348 each study for the metabolites required for certain metabolic processes. 349 We also categorised taxa by their status as potential pathogens. We matched taxon names against the can re-write F x as: (3) 360 Here, each compartment's total flux contribution (identified by a subscript: autotrophic eukaryotes = non-prokaryotic heterotrophs such as fungi or insects). We do not use the Sharpe-Schoolfield model here 366 because it does not apply to long-term thermal responses (Fig. 1), whilst for short-term responses most 367 warming as well as temperature fluctuations are expected to occur within an "operational temperature 368 range", which excludes temperatures greater than T pk (the heat-stress region) 58 . We do include any 369 potential contribution of autotrophic prokaryotes (such as cyanobacteria), as these are not expected to 370 provide a significant flux contribution to a typical terrestrial ecosystem. 371 We then use eqn. 3 to calculate the percent change in ecosystem flux due to differences in activation 372 energies of the three compartments (E ae , E hp , and E he ): where F x,2 and F x,1 are the warming-induced flux changes in ecosystems with and without differences in 374 activation energies of the compartments, respectively (the value of the heat map in Fig. 4) We calculated the emergent E of the F x,2 ecosystems (flux response to warming when prokaryotic and 383 eukaryotic thermal sensitivities differ), which is the the average of activation energies for each ecosystem 384 compartment weighted by its biomass proportion: We also calculated the emergent Q 10 of the F x,2 ecosystems, as it is a widely used measure in climate change models of carbon flux 33,59 : We chose a warming magnitude x = 10 • C for short-term responses because this at the upper end (e.g., 386 generally, at higher latitudes) of the range of daily (over 24 hrs) fluctuations that organisms experience 35 .

387
For long-term warming scenarios, we used x = 4 • C, the approximate upper end of the range for the year 388 2100 projected by the IPCC 60 .

389
The biomass proportions δ and β were varied to capture the effect of different ecosystem compositions.

390
In a typical forest ecosystem, the contribution of autotrophic to heterotrophic (mostly soil) respiration 391 has been estimated to be approximately 50% each 36 . This heterotrophic component would be comprised