Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate

Developing a thorough understanding of how ectotherm physiology adapts to different thermal environments is of crucial importance, especially in the face of global climate change. A key aspect of an organism’s thermal performance curve (TPC)—the relationship between fitness-related trait performance and temperature—is its thermal sensitivity, i.e., the rate at which trait values increase with temperature within its typically experienced thermal range. For a given trait, the distribution of thermal sensitivities across species, often quantified as “activation energy” values, is typically right-skewed. Currently, the mechanisms that generate this distribution are unclear, with considerable debate about the role of thermodynamic constraints versus adaptive evolution. Here, using a phylogenetic comparative approach, we study the evolution of the thermal sensitivity of population growth rate across phytoplankton (Cyanobacteria and eukaryotic microalgae) and prokaryotes (bacteria and archaea), 2 microbial groups that play a major role in the global carbon cycle. We find that thermal sensitivity across these groups is moderately phylogenetically heritable, and that its distribution is shaped by repeated evolutionary convergence throughout its parameter space. More precisely, we detect bursts of adaptive evolution in thermal sensitivity, increasing the amount of overlap among its distributions in different clades. We obtain qualitatively similar results from evolutionary analyses of the thermal sensitivities of 2 physiological rates underlying growth rate: net photosynthesis and respiration of plants. Furthermore, we find that these episodes of evolutionary convergence are consistent with 2 opposing forces: decrease in thermal sensitivity due to environmental fluctuations and increase due to adaptation to stable environments. Overall, our results indicate that adaptation can lead to large and relatively rapid shifts in thermal sensitivity, especially in microbes for which rapid evolution can occur at short timescales. Thus, more attention needs to be paid to elucidating the implications of rapid evolution in organismal thermal sensitivity for ecosystem functioning.


A Phylogeny reconstruction
The final tree produced by RAxML [1] and calibrated to relative time with DPPDiv [2] is shown in Fig The phylogeny generated in this study from which subtrees were extracted for comparative analyses. Colours indicate different phyla, whereas circles show the statistical support for each node, conditional to the topological constraints of the Open Tree of Life [3]. The phylogeny is available in NEXUS format at https://doi.org/10.6084/m9.figshare.12816140.v1.

B Dataset of thermal sensitivity estimates
The distributions of E and W op values across the four datasets for species included in the phylogeny are shown in Fig B. Fig C shows the distributions of thermal sensitivity estimates of r max across the six largest phyla of this study.  We first used MCMCglmm to estimate the phylogenetic heritabilities of our six TPC parameters by inferring their variance/covariance matrix, corrected for phylogeny. This allowed us to also extract the phenotypic correlation (r phe ) between E and W op and, thus, to understand the relationship between the two thermal sensitivity measures ( Fig D). Furthermore, the phenotypic correlation was broken down to its phylogenetically heritable component (r her ) and its residual component (r res ). The latter should be driven mostly by environmental effects. To understand why the relationship shown in Fig D arises, we numerically examined how W op is affected by changes in E, and the sensitivity of their relationship to changes in B 0 , T pk , and E D (Fig E).
The expected relationship of the operational niche width (W op ) and E in the Sharpe-Schoolfield model. (a) W op always decreases with E, provided that the other parameters (B 0 , T pk , and E D ) do not vary substantially and systematically with E. This is illustrated here with three arbitrary fixed combinations of the other three parameters; the curves remain practically the same, irrespective of substantial variation in these other parameters. (b) An example illustrating how W op always decreases with E when the other parameters are fixed (values shown in black). As E increases from 0.3 to 1 eV, W op decreases from 22.25 to 12.50 • C.
Besides MCMCglmm, we also estimated phylogenetic heritabilities using Rphylopars and BayesTraits and compared the resulting estimates with those of MCMCglmm ( Fig F).

MCMCglmm Rphylopars BayesTraits
Fig F. Comparison of phylogenetic heritability estimates of MCMCglmm, Rphylopars, and BayesTraits. As MCMCglmm and BayesTraits estimate phylogenetic heritability using a Bayesian approach, the plots show their posterior distributions. Instead, point estimates are shown for Rphylopars. The phylogenetic heritability estimates obtained with Rphylopars are greater than zero for all TPC parameters. While the mean phylogenetic heritability estimates of BayesTraits are generally lower than those of MCMCglmm and Rphylopars, the lower bound of the 95% Highest Posterior Density interval of BayesTraits is always greater than zero (the lowest value is at 2 · 10 −5 for ln(E) among phytoplankton). Furthermore, the distributions of phylogenetic heritabilities of prokaryotes obtained with BayesTraits are much narrower and closer to those of MCMCglmm, compared to those obtained for phytoplankton with the two programs. This is consistent with an increase in the "signal-to-noise" ratio from phytoplankton to prokaryotes, given that the latter dataset is larger (Fig B). In any case, the observed differences in the estimates of MCMCglmm and BayesTraits may arise from i) differences in the priors used by the two methods, ii) accounting (or not) for the uncertainty of each TPC parameter estimate, or from iii) differences in the approaches employed for the estimation of missing TPC parameter values. The raw data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.
To examine how the evolutionary rate of thermal sensitivity varies across the phylogeny, we fitted the stable model of trait evolution [4] (Fig 5 in the main text) but also the free model [5] (Fig G) and the Lévy model [6] ( Fig H).
Finally, to better understand how species explore the parameter space of E, W op , and T pk (whose phylogenetic heritability is ≈ 1; see Figs 2 and F), we combined our two r max datasets and divided the distributions of E, W op , and T pk into four discrete states (Fig I). Boundaries for these states were selected using the Jenks natural breaks clustering algorithm [7], as implemented in the BAMMtools R package (v. 2.1.6) [8]. To estimate the transition rates among states, we fitted the "all-rates-different" variant of the Mk model [9] with the fitMk function of the phytools R package (v. 0.6-60) [10]. Transitions in the discretized parameter space of E, W op , and T pk . The width of the edges represents the natural logarithm of the transition rate between states. Transitions between non-neighbouring states are very common for E (which captures the rise of the TPC), rare for W op (which captures both the rise and the peak of the TPC), and never observed for T pk (which captures the peak of the TPC). It is worth pointing out that T pk also exhibits the lowest transition rates between neighbouring states among the three TPC parameters. These results are consistent with the phylogenetic heritability estimates shown in Fig 2 in the main text. The raw data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.

C.2 Analyses of the dataset of phytoplankton TPCs after excluding Cyanobacteria
Removing Cyanobacteria from the phytoplankton dataset led to qualitatively identical results in our phylogenetic analyses (Figs J and K). Phylogenetic heritabilities of TPC parameters of eukaryotic phytoplankton. The main differences between these results and those using the entire phytoplankton dataset (Fig 2A in the main text) were that, here, ln(E) and ln(B pk ) have slightly higher/lower phylogenetic heritabilities respectively. The former is expected as the thermal sensitivity distribution of Cyanobacteria is very similar to that of Dinophyta (Fig 4 in the main text), despite the long evolutionary distance between them. Therefore, the exclusion of Cyanobacteria would necessarily increase the phylogenetic heritability of thermal sensitivity. Similarly, ln(B pk ) in prokaryotes is more phylogenetically heritable than in phytoplankton (Fig 2 in the main text), explaining the further decrease in its phylogenetic heritability when Cyanobacteria are excluded. The data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.

C.3 Analyses of the net photosynthesis rate and respiration rate TPC datasets
The visualization of the evolution of thermal sensitivity of net photosynthesis rate and respiration rate revealed similar patterns to those of the thermal sensitivity of r max (Fig 6 in the main text). Thermal sensitivity values do not evolve gradually and tightly around a central value ( θ), but explore large parts of the parameter space due to bursts of rapid evolution.

D.2 Fitted models using latitude as a continuous predictor
We rejected models that had one or more non-intercept coefficients with a 95% HPD interval that included zero. We then used DIC to identify the most appropriate model among those remaining. Distribution of minimum generation time estimates for phytoplankton and prokaryotes. Data points were obtained by taking the inverse of all B pk estimates. Horizontal bars represent the phylogenetically-corrected median values, i.e., the inverse of the intercept of B pk from the multi-response regression models that we fitted with MCMCglmm (see the "Estimation of phylogenetic heritability for all TPC parameters using MCMCglmm, Rphylopars, and BayesTraits" subsection of the Methods in the main text). The data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.
F List of nucleotide sequences used for phylogeny reconstruction Table D. Species names and Accession IDs of small subunit rRNA gene sequences that were used in this study.  Table E. Species names and Accession IDs of cbbL/rbcL gene sequences that were used in this study.