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.

Introduction

Current climate change projections suggest that the average global temperature in 2100 will be higher than the average of 1986–2005 by 0.3 °C–4.8 °C [1], coupled with an increase in temperature fluctuations in certain areas [2]. Therefore, it is now more important than ever to understand how temperature changes affect biological systems, from individuals to whole ecosystems. At the level of individual organisms, temperature affects functional traits in the form of the thermal performance curve (TPC). Typically, this TPC, especially when the trait is a rate (e.g., respiration rate, photosynthesis, growth), takes the shape of a negatively skewed unimodal curve (Fig 1) [3, 4]. The curve increases (approximately) exponentially to a maximum (Tpk) and then also decreases exponentially, with the fall being steeper than the rise. Understanding how various aspects of the shape of this TPC adapt to a changing thermal environment is crucial for predicting how rapidly organisms can respond to climate change.

Fig 1

The TPC of ectotherm metabolic traits, as described by the Sharpe-Schoolfield model [5].

(A) Tpk (K) is the temperature at which the curve peaks, reaching a maximum height that is equal to Bpk (in units of trait performance). E and E (eV) control how steeply the TPC rises and falls, respectively. B0 (in units of trait performance) is the trait performance normalised at a reference temperature (Tref) below the peak. In addition, Wop (K), the operational niche width of the TPC, can also be calculated a posteriori as the difference between Tpk and the temperature at the rise of the TPC where B(T) = 0.5 ⋅ Bpk. We note that we use Wop instead of a metric that captures the entire TPC width because previous studies have shown that species generally experience temperatures below Tpk [6, 7]. Thus, Wop is a measure of the thermal sensitivity of the trait near typically experienced temperatures. (B) TPCs of individual- and population-level traits (such as rmax) are usually well described by the Sharpe-Schoolfield model. The raw data for panel B are available at https://doi.org/10.6084/m9.figshare.12816140.v1. TPC, thermal performance curve.

According to the Metabolic Theory of Ecology (MTE) as well as a large body of physiological research, the shape of the TPC is expected to reflect the effects of temperature on the kinetics of a single rate-limiting enzyme involved in key metabolic reactions [5, 8–11]. Under this assumption, the rise in trait values up to Tpk can be mechanistically described using the Boltzmann-Arrhenius equation:

Here, B is the value of a biological trait, B0 is a normalisation constant—that includes the effect of cell or body size—which gives the trait value at a reference temperature (Tref), T is temperature (in K), k is the Boltzmann constant (8.617⋅10−5 eV ⋅ K−1), and E (eV) is the thermal sensitivity of the trait at the rising component of the TPC up to Tpk. Because Tpk tends to be higher than the mean environmental temperature [6, 7, 12], E represents the thermal sensitivity within the organism’s typically experienced thermal range.

Early MTE studies argued that, because of strong thermodynamic constraints, adaptation will predominantly result in changes in B0, whereas E will remain almost constant across traits (e.g., respiration rate, rmax), species, and environments around a range of 0.6–0.7 eV [8-10]. The latter assumption is referred to in the literature as universal temperature dependence (UTD). This restricted range of values that E can take is centered on the putative mean activation energy of respiration (≈ 0.65 eV). A notable exception to the UTD is photosynthesis rate, which is expected to have a lower E value of ≈ 0.32 eV, reflecting the activation energy of the rate-limiting steps of photosynthesis [13].

The existence of a UTD has been strongly debated. From a theoretical standpoint, critics of the UTD have argued that the Boltzmann-Arrhenius model is too simple to mechanistically describe the complex physiological mechanisms of diverse organisms [3, 14–16] and is inadequate for describing TPCs emerging from the interaction of multiple factors, and not just the effects of temperature on enzyme kinetics. That is, the E calculated by fitting the Boltzmann-Arrhenius model to biological traits is an emergent property that does not directly reflect the activation energy of a single rate-limiting enzyme. For example, a fixed thermal sensitivity for net photosynthesis rate is not realistic because it depends on the rate of gross photosynthesis as well as photorespiration, which is in turn determined not only by temperature but also by the availability of CO2 in relation to O2 [17].

Indeed, there is now overwhelming empirical evidence for variation in E (thermal sensitivity) far exceeding the narrow range of 0.6–0.7 eV, with such variation being, to an extent, taxonomically structured [12, 18–23]. Furthermore, the distribution of E values across species is typically not Gaussian but right-skewed. If we assume that E is nearly constant across species—and therefore that variation in E is mainly due to measurement error—such skewness could be the outcome of the proximity of the E distribution to its lower boundary (0 eV). In that case, however, we would expect a high density of E values close to 0 eV, but such a pattern has not been observed [18]. Both the deviations from the MTE expectation of a heavily restricted range for E and the shape of its distribution have been argued to be partly driven by adaptation to local environmental factors by multiple studies. These include selection on prey to have lower thermal sensitivity than predators (the “thermal life-dinner principle”) [18], adaptation to temperature fluctuations within and/or across generations [3, 21, 24–26], and adaptive increases in carbon allocation or use efficiency due to warming [27-30].

In general then, adaptive changes in the TPCs of underlying (fitness-related) traits are expected to influence the TPCs of higher-order traits such as rmax, resulting in deviations from a UTD. Therefore, understanding how the thermal sensitivity of rmax and its distribution evolves is particularly important, as it may also yield useful insights about the evolution of the TPCs of underlying physiological traits (e.g., respiration rate, photosynthesis rate, and carbon allocation efficiency). Indeed, systematic shifts in the thermal sensitivity of fundamental physiological traits have been documented [27, 31–33], albeit not through comparative analyses of large datasets.

In particular, phylogenetic heritability—the extent to which closely related species have more similar trait values than species chosen at random—can provide key insights regarding the evolution of thermal sensitivity. A phylogenetic heritability of 1 indicates that the evolution of the trait across the phylogenetic tree is indistinguishable from a random walk (Brownian motion) in the parameter space. Note that this does not necessarily indicate that the trait evolves neutrally, as it may be under selection towards a nonstationary optimum that itself performs a random walk [34]. In contrast, a phylogenetic heritability of 0 indicates that trait values are independent of the phylogeny. This is the case either because (i) the trait is practically invariant across species and any variation is due to measurement error, or (ii) the evolution of the trait is very fast and with frequent convergence (i.e., independent evolution of similar trait values by different lineages). It is worth clarifying that rapid trait evolution that does not result in convergence (e.g., when major clades are extremely separated in the parameter space) will not lead to a complete absence of phylogenetic heritability. Phylogenetic heritabilities between 0 and 1 reflect deviations from Brownian motion (e.g., due to occasional patterns of evolutionary convergence). Among phytoplankton, measures of thermal sensitivity of rmax (E and Wop) have previously been shown to exhibit intermediate phylogenetic heritability [35]. This indicates that, among phytoplankton, thermal sensitivity is not constant but evolves along the phylogeny, albeit not as a purely random walk in trait space, reflecting either thermodynamically constrained evolution or rapid evolution in response to selection.

To understand (i) how variation in thermal sensitivity accumulates across multiple autotroph and heterotroph groups and (ii) whether its distribution is shaped by environmental selection, here we conduct a thorough investigation of the evolutionary patterns of thermal sensitivity, focusing particularly on rmax. Using a phylogenetic comparative approach, we test the following hypotheses:

1. Thermal sensitivity does not evolve across species and any variation is noise-like

In this scenario, thermodynamic constraints would force E to be distributed around a mean of 0.65 eV (or 0.32 eV in the case of photosynthesis), with deviations from the mean being mostly due to measurement error. Depending on the magnitude of the error, the E distribution would either be approximately Gaussian (little measurement error) or non-Gaussian with a high density near 0 eV (substantial measurement error). This hypothesis agrees with the UTD concept of early MTE studies. If this hypothesis holds, thermal sensitivity would have 0 phylogenetic heritability and would not vary systematically across different environments.

2. Thermal sensitivity evolves gradually across species but tends to revert to a key central value, without ever moving very far from it

This hypothesis is also consistent with the UTD assumption, as it is a relaxed variant of hypothesis 1. Here, small deviations from the central tendency of 0.65 eV (or 0.32 eV) are possible, as they would reflect adaptation of species’ enzymes to certain ecological lifestyles or niches. Therefore, thermal sensitivity would be weakly phylogenetically heritable. Thermodynamic constraints would prevent large deviations from the central tendency.

3. Thermal sensitivity evolves in other ways

This is an “umbrella” hypothesis that encompasses multiple subhypotheses that do not invoke the UTD assumption. For example, a key central tendency (thermodynamic constraint) may still exist, but its influence would be very weak, allowing for a wide exploration of the parameter space away from it. In this case, changes in thermal sensitivity could be the outcome of adaptation to different thermal environments. Another subhypothesis is that clades differ systematically in the rate at which thermal sensitivity evolves, due to the occasional emergence of evolutionary innovations. Thus, clades with high evolutionary rates would be able to better explore the parameter space of thermal sensitivity (i.e., through large changes in E and Wop values), compared to low-rate clades in which thermal sensitivity would evolve more gradually. A third subhypothesis is that evolution may favour individuals (and metabolic variants) that are relatively insensitive to temperature fluctuations. In that case, the central tendency of E would not be stationary but moving towards lower values with evolutionary time. It is worth clarifying that these 3 subhypotheses are not necessarily mutually exclusive.

Results

Dataset sources

We combined 2 preexisting datasets of rmax TPCs, spanning 380 phytoplankton species (a polyphyletic group that includes prokaryotic Cyanobacteria and eukaryotic phyla such as Dinophyta) [35] and 272 prokaryote species (bacteria and archaea) [32]. In addition, we also collected 2 TPC datasets of traits that underlie rmax: net photosynthesis and respiration rates of algae and aquatic and terrestrial plants (221 and 201 species, respectively) [30]. We used these 2 smaller datasets to understand whether the evolutionary patterns of thermal sensitivity differ between (i) higher-order traits and (ii) traits that are more tightly linked to organismal physiology. Trait values were typically measured under nutrient-, light-, and CO2-saturated conditions (where applicable), after acclimation to each experimental temperature.

To investigate the evolution of measures of thermal sensitivity across species, we reconstructed the phylogeny of as many species in the 4 datasets as possible, from publicly available nucleotide sequences of (i) the small subunit rRNA gene from all species groups and the (ii) cbbL/rbcL gene from photosynthetic prokaryotes, algae, and plants (see the Methods section). We managed to obtain small subunit rRNA gene sequences from 537 species and cbbL/rbcL sequences from 208 of them (Tables D and E in S1 Appendix).

TPC parameters were quantified for each species/strain present in the phylogeny using the Sharpe-Schoolfield model (see Fig 1 and the Methods section). The resulting estimates of E (the slope of the rise of the TPC) and Wop (the operational niche width of the TPC) were found to be right-skewed (Fig B in S1 Appendix) as has been shown previously [18, 21]. Furthermore, we did not detect a disproportionately high density of thermal sensitivity values near the lower boundary of E (0 eV), as we would expect if all variation was due to strong measurement error around a true value of, e.g., 0.65 eV. Thus, these results are not consistent with the hypothesis of a nearly invariant thermal sensitivity (hypothesis 1).

Phylogenetic comparative analyses

We next investigated the evolutionary patterns of thermal sensitivity. Given that the main focus of this study was to investigate how the thermal sensitivity of rmax (a direct measure of fitness) evolves, most of the following comparative analyses were performed on our 2 large TPC datasets (rmax of phytoplankton and prokaryotes). Besides this, the sample sizes of the 2 smaller datasets would be inadequate for obtaining robust results for many of our analyses. If an analysis makes use of all 4 datasets, this is explicitly stated.

An issue that is worth mentioning is the overlap between the datasets of phytoplankton and prokaryotic TPCs, given that both of them include Cyanobacteria. To address this, we kept Cyanobacteria as part of the phytoplankton dataset (due to their functional similarity) and did not include them in analyses of prokaryotes. We also examined whether our results were mainly driven by the long evolutionary distance between Cyanobacteria and eukaryotic phytoplankton by repeating all phytoplankton analyses after removing Cyanobacteria (see subsection C.2 in S1 Appendix).

Estimation of phylogenetic heritability

As TPC parameters capture different features of the shape of the same curve, it is likely that some of them may covary [35]. To account for this in the estimation of phylogenetic heritability, we fitted a multiresponse phylogenetic regression model using the MCMCglmm R package (version 2.26) [36] in which all TPC parameters formed a combined response. To compare the phylogenetic heritabilities of TPC parameters between planktonic photosynthetic autotrophs and other microbes (autotrophs and heterotrophs), we fitted the model separately to our 2 large TPC datasets: rmax of phytoplankton and prokaryotes. To satisfy the assumption of models of trait evolution that the change in trait values is normally distributed, we transformed all TPC parameters so that their distributions would be approximately Gaussian (see Fig 2). To ensure that the resulting phylogenetic heritability estimates did not merely reflect the priors that were used in the MCMCglmm analysis, we also estimated phylogenetic heritabilities using the R package Rphylopars (version 0.2.12) [37] and BayesTraits (version 3.0.2) [38]. The main difference between these 2 and MCMCglmm is that Rphylopars and BayesTraits cannot account for the covariance among TPC parameters.

Fig 2

Moderate to strong phylogenetic heritability can be detected in all TPC parameters, across phytoplankton and prokaryotes.

The 3 circles of each radar chart correspond to phylogenetic heritabilities of 0, 0.5, and 1. Mean phylogenetic heritability estimates—as inferred with MCMCglmm—are shown in purple, whereas the 95% HPD intervals are in dark grey. Note that we transformed all TPC parameters so that their statistical distributions would be approximately Gaussian. The data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1. HPD, highest posterior density; TPC, thermal performance curve.

The MCMCglmm analysis revealed the presence of non-negligible phylogenetic heritability in measures of thermal sensitivity (E and Wop), as well as all other TPC parameters, across phytoplankton (including or excluding Cyanobacteria) and prokaryotes (Fig 2 and Fig J in S1 Appendix). In particular, the phylogenetic heritability estimates of ln(E) and ln(Wop) were statistically different from both 0 and 1, indicating that the 2 TPC parameters evolve across the phylogeny but not in a purely random (Brownian motion) manner. It is worth stressing that even the lower bounds of the 95% highest posterior density (HPD) intervals of ln(E) and ln(Wop) were far greater than 0, allowing us to completely rule out the possibility that all variation in thermal sensitivity is due to measurement error. In general, TPC parameters exhibit a similar phylogenetic heritability between the 2 species groups. The only major exception is ln(Bpk), which is considerably more heritable among prokaryotes than among phytoplankton. This difference in phylogenetic heritability most likely reflects the strength of the positive correlation between Bpk and Tpk (a “hotter is better” pattern) in the 2 groups. More precisely, Tpk, which has a phylogenetic heritability of ≈1, is more strongly correlated with Bpk among prokaryotes [32] than among phytoplankton [35], possibly due to the differences in their cellular physiology. For example, phytoplankton growth rate depends on the interplay among the processes of photosynthesis, respiration, and cell maintenance, whose thermal sensitivities can strongly differ [30]. Qualitatively similar results were obtained from the estimation of phylogenetic heritabilities with Rphylopars and BayesTraits (S1 Appendix, Fig F). Overall, these results serve as further evidence that hypothesis 1 (that thermal sensitivity does not vary across species) can clearly be rejected.

Partitioning of thermal sensitivity across the phylogeny

To understand why thermal sensitivity has a low to intermediate phylogenetic heritability, we examined how clades throughout the phylogeny explore the parameter space (of and op) using a disparity-through-time analysis [39, 40]. At each branching point of the phylogeny, mean subclade disparity is calculated as the average squared Euclidean distance among trait values within the subclades, divided by the disparity of trait values across the entire tree. Mean subclade disparity values close to 0 indicate that the mean of the trait variances within subclades is much lower than the variance of trait values across the entire phylogeny. When the opposite occurs, the mean subclade disparity will be close to 1 or even higher. The resulting disparity line is then compared to the null expectation, i.e., an envelope of disparities obtained from simulations of Brownian motion on the same tree. Through the comparison of the observed trait disparity with the null expectation, it is possible to identify the periods of evolutionary time during which mean subclade disparity is higher or lower than expected under Brownian motion. Higher than expected subclade disparity indicates that clades converge in trait space, whereas lower than expected subclade disparity indicates that clades occupy distinct areas of parameter space. The latter pattern is consistent with an adaptive radiation, in which an initial period of rapid trait evolution is typically followed by a deceleration of the evolutionary rate as ecological niches become filled [41, 42]. Frequent episodes of higher than expected subclade disparity (evolutionary convergence) in thermal sensitivity or segregation of major clades in the parameter space would be consistent with hypothesis 3.

The mean subclade disparity of thermal sensitivity measures was considerably higher than expected near the present, highlighting an increasing overlap in the parameter space of thermal sensitivity among distinct clades (Fig 3 and Fig K in S1 Appendix). This pattern of increasing clade-wide convergence in thermal sensitivity is also apparent when comparing the thermal sensitivity distributions of different phyla (Fig 4 and Fig C in S1 Appendix). For example, the distributions of E and Wop of Proteobacteria and Bacillariophyta have similar shapes and central tendencies, despite the long evolutionary distance that separates the 2 phyla. This high convergence in thermal sensitivity space by diverse lineages suggests that variation in the 2 TPC parameters is mainly driven by adaptation to local environmental conditions, irrespective of species’ evolutionary history. In other words, it is likely that particular thermal strategies (e.g., having low thermal sensitivity) may yield significant fitness gains in certain environments (e.g., those with strong temperature fluctuations that occur predominantly across—rather than within—generations [24, 25]), leading to convergent evolution of thermal sensitivity. It is worth noting that these disparity patterns are not an artefact of a potentially inaccurate tree topology, as higher than expected subclade disparity occurs mainly near the present, where tree nodes have generally high statistical support (S1 Appendix, Fig A).

Fig 3

Change in mean subclade disparity in thermal sensitivity through time.

Shaded regions represent the 95% confidence interval of the resulting trait disparity from 10,000 simulations of random Brownian evolution on each respective subtree (subset of the entire phylogeny). The dashed line stands for the median disparity across simulations, whereas the solid line is the observed trait disparity. The latter is plotted from the root of the tree (t = 0) until the most recent internal node. The reported P values were obtained from the rank envelope test [40], whose null hypothesis is that the trait follows a random walk in the parameter space. Note that instead of a single value, a range of P values is produced for each panel, due to the existence of ties. In general, species from evolutionarily remote clades tend to increasingly overlap in thermal sensitivity space (mean subclade disparity exceeds that expected under Brownian motion) with time. The raw data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.

Fig 4

Distributions of thermal sensitivity estimates of rmax for the largest (most species-rich) phyla of this study.

In general, more variation can be observed within than among phyla. The data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.

Mapping the evolutionary rate on the phylogeny

We next investigated whether clades systematically differ in their evolutionary rate for thermal sensitivity (part of hypothesis 3). To this end, we examined the variation in the evolutionary rate of thermal sensitivity measures across the phylogeny by fitting 3 extensions of the Brownian motion model: the free model [43], the stable model [44], and the Lévy model [45]. Under the free model, the trait takes a random walk in the parameter space (Brownian motion) but with an evolutionary rate that varies across branches. The stable model can be seen as a generalisation of the free model, as the evolutionary change in trait values is sampled from a heavy-tailed stable distribution, of which the Gaussian distribution (assumed under Brownian motion) is a special case. Thus, the stable model should provide a more accurate representation of evolutionary rate variation, as it is better able to accommodate jumps in parameter space towards rare and extreme trait values. Finally, the Lévy model represents evolution under Brownian motion combined with occasional episodes of rapid trait change.

The results were robust to the choice of model used for inferring evolutionary rates (Fig 5, Figs G and H in S1 Appendix). Rate shifts tend to occur sporadically throughout the phylogeny and especially in late-branching lineages, without being limited to particular clades. This pattern suggests that there is little systematic variation in the evolutionary rate of thermal sensitivity among clades, with sudden bursts of trait evolution arising in parallel across evolutionarily remote lineages.

Fig 5

Variation in the evolutionary rate of thermal sensitivity across the phylogeny.

Rates were estimated by fitting the stable model of trait evolution to each dataset and were then normalised between 0 and 1. Most branches exhibit relatively low rates of evolution (orange), whereas the highest rates (red and brown) are generally observed in late-branching lineages across different clades. The raw data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.

Visualization of trait evolution as a function of time, and test for directional selection

To further describe the evolution of thermal sensitivity, we visualized the E and Wop values from the root of each subtree until the present day, across all 4 TPC datasets, using the phytools R package (version 0.6–60) [46]. Ancestral states—and the uncertainty around them—were obtained from fits of the stable model of trait evolution, as described in the previous subsection. The visualization allowed us to test hypothesis 2, i.e., that thermal sensitivity evolves around a central tendency of 0.65 eV (or 0.32 eV), with large deviations from this value reverting quickly back to it. To this end, and to also test the hypothesis of directional selection towards lower thermal sensitivity (part of hypothesis 3), we used the following model:

ln(E) values (those from extant species and ancestral states inferred with the stable model) were regressed against a central value (ln()) and a slope that captures a putative linear trend towards lower/higher values with relative time, t. The same model was also fitted to ln(Wop). The regressions were performed with MCMCglmm and were corrected for phylogeny as this resulted in lower deviance information criterion (DIC) [47] values than those obtained from non-phylogenetic variants of the models. More precisely, we executed 2 MCMCglmm chains per regression for a million generations, sampling every thousand generations after the first hundred thousand.

This analysis (Fig 6) did not provide support for the hypothesis of strongly constrained adaptive evolution around a single key central value (hypothesis 2). Instead, lineages explore large parts of the parameter space, often moving rapidly towards the upper and lower bounds (i.e., 0 and 4 eV), without reverting back to the presumed central tendency (e.g., see the clade denoted by the arrow in Fig 6D). The estimated central values for E of the two rmax datasets were much higher than the MTE expectation, and, in the case of prokaryotes (Fig 6B), the 95% HPD interval did not include 0.65. Similarly, the inferred central values for E of net photosynthesis rate and respiration rate (0.52 eV and 2.06 eV, respectively; Fig La,b in S1 Appendix) were both higher than 0.32 and 0.65 eV. The slope parameter that would capture the effects of directional selection in thermal sensitivity (part of hypothesis 3) was not statistically different from 0 for any dataset.

Fig 6

Projection of the phylogeny into thermal sensitivity versus time space.

The values of ancestral nodes were estimated from fits of the stable model. Yellow lines represent the median estimates, whereas the 95% credible intervals are shown in red. is the estimated central tendency for each panel, whereas the existence of a linear trend towards lower/higher values is captured by the reported slope. Parentheses stand for the 95% HPD intervals for and the slope. All estimates were obtained for ln(E) and ln(Wop), but the parameters are shown here in linear scale. The inset figures show the density distributions of E and Wop values of extant species in the dataset. The arrow in panel D shows an example of a whole clade shifting towards high Wop values, without being attracted back to . The raw data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1. HPD, highest posterior density.

Latitudinally structured variation in thermal sensitivity

All our analyses so far converge on one conclusion: that the evolution of the thermal sensitivities of fitness-related traits can be rapid and largely independent of the evolutionary history of each lineage. This suggests that certain environments may select for particular values of thermal sensitivity. To identify environmental adaptation in thermal sensitivity, we tested for latitudinal variation in it using the combination of all 4 TPC datasets. Specifically, the increase in temperature fluctuations from low to intermediate absolute latitudes is expected to increasingly select for thermal generalists (lower E and higher Wop values) [3, 48, 49]. At high latitudes, however, temperature fluctuations may further increase or progressively decrease, depending on environment type (marine versus terrestrial) and differences between the 2 hemispheres [3, 48, 49]. In any case, the overwhelming majority of our thermal sensitivity estimates belonged to species/strains from low and intermediate latitudes (S1 Appendix, Fig M), enabling us to investigate the hypothesized gradual transition towards lower thermal sensitivity from the equator to intermediate latitudes.

Latitude indeed explained a significant amount of variation in E (which declined as expected) but not in Wop (Fig 7 and Fig N in S1 Appendix, Tables A and B in S1 Appendix). The E estimates of rmax, net photosynthesis rate, and respiration rate differed statistically in their intercepts but not in their slopes against latitude, although the latter could be an artefact of the small sample size. This result suggests that latitude could influence the E values of not only rmax but also other traits across various species groups. Dividing latitude into 3 bins (i.e., low, intermediate, and high absolute latitudes) and comparing their E distributions yielded similar conclusions (S1 Appendix, Fig O, Table C).

Fig 7

E values weakly decrease with absolute latitude.

23% of the variance is explained by latitude and trait identity, which increases to 58% if species identity is added as a random effect on the intercept. Note that values on the vertical axis increase exponentially. The data underlying this figure are available at https://doi.org/10.6084/m9.figshare.12816140.v1.

We also tested for a possible latitudinal clade age bias, which could arise if certain clades originated in particular latitudes and only much later expanded to other areas [50, 51]. For this, we performed a Mantel test [52] to estimate the correlation between phylogenetic distance and latitudinal distance for the 2 largest groups of our study (phytoplankton and prokaryotes). No such bias was detected for phytoplankton (r = 0.04, P = 0.114), whereas, for prokaryotes, the correlation was statistically supported but very weak (r = 0.11, P = 0.002). This result indicates that neither species group is characterised by very strong dispersal limitation throughout its evolutionary history.

Discussion

In this study, we have performed a thorough analysis of the evolution of the thermal sensitivities of rmax in phytoplankton and prokaryotes and its 2 key underlying physiological traits in plants (net photosynthesis rate and respiration rate). To achieve this, we formulated and tested 3 alternative hypotheses that represent different views expressed in the literature regarding the impact of thermodynamic constraints on the evolution of thermal sensitivity of fitness-related traits.

The first hypothesis was that the activity of a single key rate-limiting enzyme of respiration or photosynthesis directly determines the performance of physiological traits [9, 13] and emergent proxies for fitness (such as rmax) [53, 54] (the UTD assumption). As a result, thermal sensitivity should be strictly constant across traits, species, and environments. This hypothesis was first introduced in early papers that described the MTE [8-10]. In contrast to the UTD expectation, we detected substantial variation in thermal sensitivity, within and across traits and species groups (Fig B in S1 Appendix). Furthermore, the distribution of E (slope of the rising part of the TPC) values did not exhibit an inflated density near its lower boundary (around 0 eV), as we would expect if all variation in thermal sensitivity was due to measurement error. The rejection of hypothesis 1 was additionally supported by our finding that thermal sensitivity is phylogenetically heritable across phytoplankton and prokaryotes (Fig 2).

Our second hypothesis was that thermal sensitivity evolves across species but remains close to a key value imposed by strong (but not insurmountable) thermodynamic constraints. We tested this hypothesis using a series of phylogenetic comparative analyses which revealed that the evolution of thermal sensitivity is characterised by an increasing overlap in parameter space by evolutionarily remote lineages (Figs 3 and 4) due to bursts of rapid evolution (Fig 5). Additionally, visualisation of thermal sensitivity evolution through time (Fig 6 and Fig L in S1 Appendix) showed that thermal sensitivity can rapidly move away from its presumed central value without being strongly attracted back to it (e.g., see the arrow in Fig 6D). In conclusion, these results lead us to reject hypothesis 2, i.e., that thermal sensitivity evolves under very strong thermodynamic constraints.

Our final hypothesis was that thermal sensitivity evolves in an adaptive manner and that, even if a central tendency exists, its influence on thermal sensitivity evolution is weak. This hypothesis was supported by the results of all phylogenetic comparative analyses and by our detection of a systematic relationship between E and latitude. The latter is likely driven by the increase in temperature fluctuations from the equator to intermediate latitudes and agrees with the expectation that thermally variable environments should select for phenotypes with low thermal sensitivity, and vice versa [3, 21, 24, 25]. That a similar latitudinal effect could not be detected on Wop (the operational niche width of the TPC) is possibly because of the much smaller sample size available for this, combined with the fact that Wop is nonlinearly related to E (S1 Appendix, Fig D). More precisely, in the Sharpe-Schoolfield model, Wop will necessarily decrease as E increases, provided that changes in E are not strongly associated with changes in other TPC parameters (e.g., Tpk or E; Fig E in S1 Appendix). Indeed, a previous study showed that E correlates systematically only with Wop [35]. In any case, E is arguably a more meaningful measure of thermal sensitivity than Wop because the latter assumes that species mainly experience temperatures close to Tpk, while E captures the entire rise of the TPC. Besides temperature fluctuations, a decrease in E with absolute latitude could also be explained by the metabolic cold adaptation hypothesis [55-58]. According to it, cold-adapted species should evolve lower thermal sensitivities (as well as higher B0 values; see Fig 1) to maintain sufficient trait performance at very low temperatures. As our datasets do not possess the necessary resolution (especially at high latitudes; Figs M and O in S1 Appendix) for differentiating between these 2 alternative (and non-mutually exclusive) processes, this question remains to be addressed by future research.

Overall, a set of novel mechanistic explanations of TPC evolution emerge from our comparison of phylogenetic heritabilities of TPC parameters (Fig 2 and Fig F in S1 Appendix). Contrary to E and Wop, which have low to intermediate phylogenetic heritabilities, Tpk is almost perfectly phylogenetically heritable and evolves relatively gradually (i.e., without large jumps in parameter space; see Fig I in S1 Appendix). Thus, we expect TPCs to adapt to different thermal environments through both gradual changes in Tpk and discontinuous changes in E. Gradual changes in Tpk may be achieved through evolutionary shifts in the melting temperature of enzymes, i.e., the temperature at which 50% of the enzyme population is deactivated [59, 60]. In contrast, changes in thermal sensitivity may be the outcome of (i) evolution of enzymes with different heat capacities [60-62], (ii) changes in the plasticity of cellular membranes [3, 63], or even (iii) restructuring of the underlying metabolic network [64].

Fundamental differences in the selection mechanisms underlying the evolution of Tpk and E may also explain the difference in evolutionary patterns between them. Specifically, both the mean environmental temperature (to which Tpk responds [7, 35]) and the temperature fluctuations (to which E responds [3, 21, 24–26, 35]) vary systematically from the equator to intermediate latitudes. We hypothesize that a species adapted to low temperatures is unlikely to adapt to a high-temperature environment rapidly enough (i.e., through a large increase in Tpk) as it is pushed to its thermal tolerance limits [65, 66]. In contrast, a species adapted to a fluctuating thermal environment (i.e., with a low E value) should be able to survive in more thermally stable conditions without much cost, becoming a thermal specialist (i.e., with a high E value) relatively rapidly, resulting in the observed jumps in trait space when mapped on the phylogeny (Figs 3 and 6, Figs K and L in S1 Appendix).

It is worth stressing, however, that not all types of thermal fluctuations are expected to impose selection for thermal generalists. In particular, thermal generalist variants of a given species are expected to be favoured when temperature fluctuations primarily occur across generations [24, 25, 67]. In contrast, moderate to strong thermal variation within generations would lead to selection for thermal specialists, even when intergenerational fluctuations are also present. For the microbial groups of the present study, an estimate of the minimum generation time can be calculated as the inverse of the Bpk of rmax. Across our datasets of phytoplankton and prokaryotes, the minimum generation time ranges from a few minutes to 3.5 months, with phylogenetically corrected medians of ≈ 40.5 hours for phytoplankton and ≈ 3.5 hours for prokaryotes (S1 Appendix, Fig P). Given this and because the magnitudes of annual and intra-annual (e.g., monthly) thermal fluctuations increase from the equator to intermediate latitudes [3, 48, 68], most microbes from intermediate latitudes are expected to generally experience substantial intergenerational thermal fluctuations and to a much lesser extent intragenerational fluctuations. This is indeed consistent with the observed weak decline in E at intermediate latitudes compared to the equator (Fig 7 and Fig O in S1 Appendix). Nevertheless, latitude, trait identity, and species identity account for only 58% of the variance in E, indicating that adaptive shifts in E may also be driven by other factors such as biotic interactions [18, 69, 70]. A systematic identification of drivers of thermal sensitivity as well as the magnitude of their respective influence could be the focus of future studies.

For the thermal sensitivity of rmax in particular, the observed patterns of discontinuous evolution likely reflect the evolution of TPCs of underlying physiological traits on which it depends. For example, in populations of photosynthetic cells, shifts in the thermal sensitivity of any or all of photosynthesis rate, respiration rate, and carbon allocation efficiency can induce large changes in the E of rmax [30]. Indeed, we observed large adaptive shifts in thermal sensitivity even for fundamental physiological traits such as respiration rate (S1 Appendix, Fig Lb,d), contrary to the MTE expectation of strong evolutionary conservatism [8-10]. This result is in agreement with a previous study that had identified significant adaptive variation in the TPC of the specific activity of Rubisco carboxylase [31]. It remains to be seen whether a similar lack of evolutionary conservation can be detected in key enzymes of non-photosynthetic organisms. Further research is clearly also needed on how the thermal sensitivities of different traits underlying fitness interact, and the extent to which these interactions can be modified through adaptation.

Besides biological-driven variation in thermal sensitivity, “artificial” variation may also be present, hindering the recognition of real patterns. For example, E estimates can be inaccurate if trait measurements in the rise of the TPC are limited, and span too narrow a range of temperatures [12]. To address this issue, we only kept E estimates if at least 4 trait measurements were available at the rise of each TPC. Further variation in thermal sensitivity can be introduced if trait values are measured instantaneously (without allowing sufficient time for acclimation) or under suboptimal conditions (e.g., under nutrient- or light-deficient conditions). Such treatments can lead to systematic biases in the shape of the resulting TPCs, which may strongly differ from TPCs obtained after adequate acclimation and under optimal growth conditions [27, 71–74]. To avoid such biases, the datasets that we used only included TPCs that were experimentally determined after acclimation and under optimal conditions. On the other hand, maintenance of a given strain under a fixed set of experimental conditions for hundreds of generations could also lead to adaptive changes in TPC shape, due to the emergence of novel genetic mutations, as has been previously shown [26, 27]. While the strains in our dataset were not grown over such long time periods, future studies could employ experimental evolution to measure the rate of thermal sensitivity evolution over much shorter timescales than the ones in our study.

Put together, all these results yield a compelling mechanistic explanation of how evolution shapes the distribution of E and emphasize the need to consider the ecological and evolutionary underpinnings as well as implications of variation in E, as has been pointed out in a spate of recent studies [12, 18, 20, 21, 30]. In particular, our study helps explain the reason for the right skewness in the E distributions previously identified across practically all traits and taxonomic groups [12, 18, 21]. A clear explanation for this pattern has been lacking, partly because MTE posits that E should be thermodynamically constrained and thus almost invariable across species [8-10]. Our study fills this gap in understanding by showing that the distribution of E is the outcome of frequent convergent evolution, driven by the adaptation of species from different clades to similar environmental conditions. In other words, as species encounter new environments through active or passive dispersal [75-77], they face selection for particular values of thermal sensitivity, which results in (often large) shifts in E. This process explains both the low variation in E among some species groups (Fig 4) and the shape of its distribution. More precisely, the high degree of right skewness probably reflects the fact that most environments select for thermal generalists, with high E values being less frequently advantageous. Our findings have implications for ecophysiological models, which may benefit from accounting for variation in thermal sensitivity among species or individuals. This could both yield an improved fit to empirical datasets [78] and provide a more realistic approximation of the processes being studied. Finally, the existence of adaptive variation in thermal sensitivity is likely to partly drive ecological patterns at higher scales (e.g., the response of an ecosystem to warming). How differences in thermal sensitivity among species influence ecosystem function is largely unaddressed [32, 78] but highly important for accurately predicting the impacts of climate change on diverse ecosystems.

Methods

Phylogeny reconstruction and relative time calibration

We performed sequence alignment using MAFFT (version 7.123b) [79] and its L-INS-i algorithm, and we ran Noisy (version 1.5.12) [80] with the default options to identify and remove phylogenetically uninformative homoplastic sites. For a more robust phylogenetic reconstruction, we used the results of previous phylogenetic studies by extracting the Open Tree of Life [81] topology for the species in our dataset using the rotl R package [82]. We manually examined the topology to eliminate any obvious errors. In total, 497 species were present in the tree, whereas many nodes were polytomic. To add missing species and resolve polytomies, we inferred 1,500 trees with RAxML (version 8.2.9) [83] from our concatenated sequence alignment, using the Open Tree of Life topology as a backbone constraint and the General Time-Reversible model [84] with Γ-distributed rate variation among sites [85]. This model was fitted separately to each gene partition (i.e., one partition for the alignment of the small subunit rRNA gene sequences and one partition for the alignment of cbbL/rbcL gene sequences). Out of the 1,500 resulting tree topologies, we selected the tree with the highest log-likelihood and performed bootstrapping (using the extended majority-rule criterion) [86] to evaluate the statistical support for each node.

Finally, we calibrated the resulting RAxML tree to units of relative time by running DPPDiv [87] on the alignment of the small subunit rRNA gene sequences using the uncorrelated Γ-distributed rates model [88] (S1 Appendix, Fig A). For this, we used the alignment of small subunit rRNA gene sequences only, as DPPDiv can only be run on a single gene partition. We executed 2 DPPDiv runs for 9.5 million generations, sampling from the posterior distribution every 100 generations. After discarding the first 25% of samples as burn-in, we ensured that the 2 runs had converged on statistically indistinguishable posterior distributions by examining the effective sample size and the potential scale reduction factor [89, 90] for all model parameters. More precisely, we verified that all parameters had an effective sample size above 200 and a potential scale reduction factor value below 1.1. To summarise the posterior distribution of calibrated trees into a single relative chronogram, we kept 4,750 trees per run (one tree every 1,500 generations) and calculated the median height for each node using the TreeAnnotator program [91].

Sharpe-Schoolfield model fitting

To obtain estimates of the parameters of each experimentally determined TPC, we fitted the following 4-parameter variant of the Sharpe-Schoolfield model (Fig 1) [5, 35]:

This model extends the Boltzmann-Arrhenius model (Eq 1) to capture the decline in trait performance after the TPC reaches its peak (Tpk). We followed the same approach for fitting the Sharpe-Schoolfield model as Kontopoulos and coworkers [35]. Briefly, we set Tref to 0 °C because, for B0 to be biologically meaningful (see Fig 1), it needs to be normalised at a temperature below the minimum Tpk in the study. Thus, a Tref value of 0 °C allowed us to include TPCs from species with low Tpk values in the analyses. Also, as certain specific TPC parameter combinations can mathematically lead to an overestimation of B0 compared to the true value, B(Tref) [92], we manually recalculated B(Tref) for each TPC after obtaining estimates of the 4 main parameters (B0, E, Tpk, and E). For simplicity, these recalculated B(Tref) values are referred to as B0 throughout the study. Finally, Bpk and Wop were calculated based on the estimates of the 4 main parameters.

After rejecting fits with an R2 below 0.5, there were (i) 312 fits across 118 species from the phytoplankton rmax dataset, (ii) 289 fits across 189 species from the prokaryote rmax dataset, (iii) 87 fits across 38 species from the net photosynthesis rates dataset, and (iv) 34 fits across 18 species from the respiration rates dataset. Note that some species were represented by multiple fits due to the inclusion of experimentally determined TPCs from different strains of the same species or from different geographical locations. To ensure that all TPC parameters were reliably estimated, we performed further filtering based on the following criteria: (i) B0 and E estimates were rejected if fewer than 4 experimental data points were available below Tpk. (ii) Extremely high E estimates (i.e., above 4 eV) were rejected. (iii) Wop values were retained if at least 4 data points were available below Tpk and 2 after it. (iv) Two data points below and after the peak were required for accepting the estimates of Tpk and Bpk. (v) E estimates were kept if at least 4 data points were available at temperatures greater than Tpk.

Estimation of phylogenetic heritability for all TPC parameters using MCMCglmm, Rphylopars, and BayesTraits

For MCMCglmm, the methodology that we used was also identical to that of Kontopoulos and coworkers [35]. In short, we specified a phylogenetic mixed-effects model for each of the 2 large TPC datasets. The models had a combined response with all TPC parameters transformed towards normality. The uncertainty for each estimate was obtained with the delta method [93] or via bootstrapping (for ln(Wop)) and was incorporated into the model. Missing estimates in the response variables (i.e., when not all parameter estimates could be obtained for the same TPC) were modelled according to the “Missing At Random” approach [36, 94]. Regarding fixed effects, a separate intercept was specified for each TPC parameter. Species identity was treated as a random effect on the intercepts and was corrected for phylogeny through the integration of the inverse of the phylogenetic variance-covariance matrix. For each dataset, 2 Markov chain Monte Carlo chains were run for 200 million generations, and estimates of the parameters of the model were sampled every 1,000 generations after the first 20 million generations were discarded as burn-in. Tests to ensure that the chains had converged and that the parameters were adequately sampled were done as previously described.

We also estimated Pagel’s λ [95] (which is equivalent to phylogenetic heritability [96]) for each TPC parameter using Rphylopars and BayesTraits. For the latter, we executed 2 Markov chain Monte Carlo chains for 10 million generations, kept samples from the posterior every 1,000 generations after the first million, and ensured that sufficient convergence had been reached. Nevertheless, we note that our previous approach is superior because Rphylopars and BayesTraits analyse each TPC parameter separately, and thus covariances among TPC parameters are not taken into account when estimating missing values. Furthermore, these 2 methods cannot accommodate the uncertainty for each TPC parameter estimate.

Disparity-through-time analyses

We performed disparity-through-time analyses for ln(E) and ln(Wop), using the rank envelope method [40] to generate a confidence envelope from 10,000 simulations of random evolution (Brownian motion). As it is not straightforward to incorporate multiple measurements per species with this method, we selected the ln(E) or ln(Wop) estimate of the Sharpe-Schoolfield model fit with the highest R2 value per species.

Free, stable, and Lévy model fitting

We fitted the free, stable, and Lévy models of trait evolution to estimates of ln(E) and ln(Wop), using the motmot.2.0 R package (version 1.1.2) [97, 98], the stabletraits software [44], and the levolution software [45], respectively. To obtain each fit of the stable model, we executed 4 independent Markov chain Monte Carlo chains for 30 million generations, recording posterior parameter samples every 100 generations. Samples from the first 7.5 million generations were excluded, whereas the remaining samples were examined to ensure that convergence had been achieved. For fitting the Lévy model, we used the peak-finder algorithm to estimate the value of the model’s α parameter. More precisely, we set the starting value of α to 100.5, the step size to 0.5, and the number of optimizations to 5, as suggested in levolution’s documentation. We also changed the maximum number of iterations (option “-maxIterations”) to 2,000 so that the algorithm could sufficiently converge in all cases.

Investigation of a putative relationship between latitude and ln(E) and ln(Wop)

We examined the relationship of thermal sensitivity with latitude by fitting regression models with MCMCglmm to all 4 TPC datasets combined. The response variable was ln(E) or ln(Wop), whereas possible predictor variables were (i) latitude (in radian units and using a cosine transformation, as absolute latitude in degree units, or split in 3 bins of low, intermediate, and high absolute latitude; subsections D.2 and D.3 in S1 Appendix), (ii) the trait from which thermal sensitivity estimates were obtained, and (iii) the interaction between latitude and trait identity. To properly incorporate multiple measurements from the same species (where available), we treated species identity as a random effect on the intercept. We fitted both phylogenetic and non-phylogenetic variants of all candidate models. Two chains per model were run for 5 million generations each, with samples from the posterior being captured every thousand generations. We verified that each pair of chains had sufficiently converged, after discarding samples from the first 500,000 generations. To identify the most appropriate model, we first rejected models that had a nonintercept coefficient with a 95% HPD interval that included 0. We then selected the model with the lowest mean DIC value. To report the proportions of variance explained by the fixed effects (Varfixed), by the random effect (Varrandom), or left unexplained (Varresid), we calculated the marginal and conditional coefficients of determination [99]:

Mantel test between phylogenetic and latitudinal distance matrices

We used the R package ade4 (version 1.7–13) [100] to infer the correlation of phylogenetic distance with latitudinal distance across phytoplankton and prokaryotes using the Mantel test. To generate the P values, we set the number of permutations to 9,999.

Supplementary material.

(PDF)

Click here for additional data file.

9 Jan 2020

Dear Dr Kontopoulos,

Thank you for submitting your manuscript entitled "Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate" for consideration as a Research Article by PLOS Biology.

Your manuscript has now been evaluated by the PLOS Biology editorial staff and I am writing to let you know that we would like to send your submission out for external peer review.

However, before we can send your manuscript to reviewers, we need you to complete your submission by providing the metadata that is required for full assessment. To this end, please login to Editorial Manager where you will find the paper in the 'Submissions Needing Revisions' folder on your homepage. Please click 'Revise Submission' from the Action Links and complete all additional questions in the submission questionnaire.

Please re-submit your manuscript within two working days, i.e. by Jan 11 2020 11:59PM.

Login to Editorial Manager here: https://www.editorialmanager.com/pbiology

During resubmission, you will be invited to opt-in to posting your pre-review manuscript as a bioRxiv preprint. Visit http://journals.plos.org/plosbiology/s/preprints for full details. If you consent to posting your current manuscript as a preprint, please upload a single Preprint PDF when you re-submit.

Once your full submission is complete, your paper will undergo a series of checks in preparation for peer review. Once your manuscript has passed all checks it will be sent out for review.

Feel free to email us at plosbiology@plos.org if you have any queries relating to your submission.

Kind regards,

Lauren A Richardson, Ph.D

Senior Editor

PLOS Biology

4 Feb 2020

Dear Dr Kontopoulos,

Thank you very much for submitting your manuscript "Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate" for consideration as a Research Article at PLOS Biology. Your manuscript has been evaluated by the PLOS Biology editors, an Academic Editor with relevant expertise, and by several independent reviewers.

As you will read, the reviewers appreciated many aspects of your study. However, Reviewer #3 raises a number of key issues with the phylogenetic methods employed. Further analyses are needed in a revision to support these conclusions.

In light of the reviews (below), we will not be able to accept the current version of the manuscript, but we would welcome re-submission of a much-revised version that takes into account the reviewers' comments. We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent for further evaluation by the reviewers.

We expect to receive your revised manuscript within 2 months.

Please email us (plosbiology@plos.org) if you have any questions or concerns, or would like to request an extension. At this stage, your manuscript remains formally under active consideration at our journal; please notify us by email if you do not intend to submit a revision so that we may end consideration of the manuscript at PLOS Biology.

**IMPORTANT - SUBMITTING YOUR REVISION**

Your revisions should address the specific points made by each reviewer. Please submit the following files along with your revised manuscript:

1. A 'Response to Reviewers' file - this should detail your responses to the editorial requests, present a point-by-point response to all of the reviewers' comments, and indicate the changes made to the manuscript.

*NOTE: In your point by point response to the reviewers, please provide the full context of each review. Do not selectively quote paragraphs or sentences to reply to. The entire set of reviewer comments should be present in full and each specific point should be responded to individually, point by point.

You should also cite any additional relevant literature that has been published since the original submission and mention any additional citations in your response.

2. In addition to a clean copy of the manuscript, please also upload a 'track-changes' version of your manuscript that specifies the edits made. This should be uploaded as a "Related" file type.

*Re-submission Checklist*

When you are ready to resubmit your revised manuscript, please refer to this re-submission checklist: https://plos.io/Biology_Checklist

To submit a revised version of your manuscript, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' where you will find your submission record.

Please make sure to read the following important policies and guidelines while preparing your revision:

*Published Peer Review*

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*PLOS Data Policy*

Please note that as a condition of publication PLOS' data policy (http://journals.plos.org/plosbiology/s/data-availability) requires that you make available all data used to draw the conclusions arrived at in your manuscript. If you have not already done so, you must include any data used in your manuscript either in appropriate repositories, within the body of the manuscript, or as supporting information (N.B. this includes any numerical values that were used to generate graphs, histograms etc.). For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5

*Blot and Gel Data Policy*

We require the original, uncropped and minimally adjusted images supporting all blot and gel results reported in an article's figures or Supporting Information files. We will require these files before a manuscript can be accepted so please prepare them now, if you have not already uploaded them. Please carefully read our guidelines for how to prepare and upload this data: https://journals.plos.org/plosbiology/s/figures#loc-blot-and-gel-reporting-requirements

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosbiology/s/submission-guidelines#loc-materials-and-methods

Thank you again for your submission to our journal. We hope that our editorial process has been constructive thus far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Lauren A Richardson, Ph.D

Senior Editor

PLOS Biology

*****************************************************

REVIEWS:

Reviewer #1:

The research presented here aims to answer two questions: How variation in thermal sensitive accumulates across multiple autotroph and heterotroph groups (meaning, over time?), and then, whether variation in thermal sensitivity is shaped by environmental selection. This is an important question and of great interest to scientists across disciplines of biology, from physiology to ecology, and the results here are certainly novel and exciting. Variation among thermal sensitivities among taxa has certainly been demonstrated before, and so has variation in temperature dependence of metabolic rates among prokaryote and eukaryote taxonomic groups (e.g., White et al Proc Roy Soc 2012 279: 1742, Galmes et al, Photosynthesis research 2015: 123(2) 183-201.) But, a robust analysis explicitly testing hypotheses from MTE about the universality of activation energy parameters that quantitatively considers phylogenies lacking, to my knowledge. Therefore, I believe this study constitutes a substantial and important contribution to the literature. My belief would be confirmed if the authors could explicitly relate their findings to the two studies noted above (and others if I've missed them) that reach similar findings by different methods (White et al) or different findings (Galmes et al) than reported in the present study.

I should note that I do not have expertise in phylogenetic methods or analyses, and so I cannot evaluate the authors' phylogenetic approach, assumptions or interpretations here. These methods are central to their study and inferences, so I feel uncomfortable providing a strong recommendation on publication or not. I hope there are other reviewers with expertise in phylogeny to provide the editor with a full assessment of this paper.

My main comments are meant to be constructive to the authors in revising their paper to be more clear to readers not expert in this topic. As I mentioned already, I think the findings of this paper are of general interest, but the potential limitations or assumptions in the work are not clear to a non-specialist reader as the manuscript is currently written. I recommend making much more clear and explicit the set of hypotheses tested, and the assumptions inherent in the hypothesis as well as in the methods used to test them. I recommend revising the introduction to explicitly state a set of testable (rejectable) hypotheses. I find myself trying to figure out what the hypotheses are, but really it is the job of the authors to tell me. My reading of the introduction leads me to think that the hypotheses being tested here are:

H0: there is no phylogenetic signal to thermal sensitivity traits (as defined by the empirically fit sharpe-schoolfield equation). This would be consistent with MTE's UTD in its simplest form. (The authors refer to this as an MTE assumption - line 89.)

One prediction associated with this hypothesis that the authors have given is that measurement error could explain any patterns in the data, rather than phylogenetic signal. The introduction nicely explains what patterns would support or refute this prediction. I would like to see the alternative predictions introduced as clearly.

Another hypothesis I draw out of the introduction is:

H1: there is phylogenetic signal to thermal sensitivity traits, reflecting environmental conditions experienced by species or groups during their evolution. The authors could be more specific about what these patterns might be expected to be.

On the geographic variation hypotheses (H2, perhaps), there are alternate explanations that should be considered (especially because thermal variation (annual?) as a predictor was not explicitly tested or even shown. Furthermore, rationale should be given for how (annual?) thermal variation is expected to matter for organisms with such high rates of population growth and presumably evolution. Additionally, many plankton have resting stages in which their populations do not grow during harsh conditions. So the link between thermal variation and diversity of thermal traits in unicellular algae needs to be explained a bit more, and alternative need to be explored. Many ecological and evolutionary patterns and processes vary with latitude, other than thermal variation. Most notably, biodiversity and evolutionary history (tropical groups are generally older). These should be considered or at least mentioned. When this comes up in the discussion (line 279), one alternative explanation is considered (MCA). The finding that E varies with latitude but not W may (or may not) be consisted with explanations other than simply different sample sizes.

The first paragraph of the discussion references two hypotheses, and it would be stronger to see these set up more obviously in the introduction. It is there, but the writing in that first paragraph of the discussion is much more clear in terms of presenting the hypotheses tested.

Lines 184-186:

"This high convergence in thermal sensitivity space by diverse lineages suggests that variation in the two TPC parameters is mainly driven by adaptation to local environmental conditions, irrespective of the species evolutionary history". Here is a place where the strength of inference and findings would be more clear if a (full?) set of alternate hypotheses were laid out clearly. From what is written, I infer that the possible interpretations for patterns of variation in TPC parameters across species are:

- Evolutionary history (presumably independent of the environmental influence?)

- Environmental conditions (independent of evolutionary history, as noted by the authors)

- Measurement error (reject above)

- What are other possibilities - constraints on possible traits associated with physical processes (rather than genetic constraints associated with evolution)? Can alternatives be rejected? Is an interaction between evolutionary history and environment possible, and if so, how would that be diagnosed?

Line 204 - 206:

This seems like an interesting result but needs more explanation and setup. Again, a set of hypotheses or predictions might help here. As a reader, I was unprepared for the possible outcomes of this test, and therefore it was unclear what were the implications of their actual finding. It would also help if, when introducing possible outcomes, the authors could highlight key assumptions of their methods and whether those assumptions might lend uncertainty to the conclusions. I think this would help a general biologist reader of PLoS biology who does not have expertise in phylogenetics (like myself).

The results and analysis in Figure 6 are hard to understand. It is not clear why variation in thermal sensitivity among taxa increases with time, nor how one can confidently reconstruct historical thermal sensitivities having previously concluded that they are not explained by evolutionary history but rather by environmental conditions. Does the reconstruction in figure 6 not assume that these traits are affected by evolution, and therefore contradict the previous result?

Methods:

Though I can't evaluate the phylogenetic methods, I do think the authors' general statistical approach of ensuring that their datasets do not include overlapping species, and the consideration of possible covariation among TPC parameters in their analysis are appropriate.

How much confidence do the authors have in the phylogenetic estimates of time (their relative time axis)? These uncertainties should be explained. Not being an expert, I assume it is likely that evolutionary relationships are better resolved for more recent divisions compared to those in deeper time; could the uncertainty in phylogenies associated with time explain part of the results in figure 3?

Minor comments:

There seems to be increasing evidence that the E value of photosynthesis may not universally be 0.32, as was suggested by Allen et al 2005. There has been support for this value, mostly in algae, but less in terrestrial plants. You might consider removing this sentence (line 29-31). However, this is very much unresolved. The work I've seen that convinced me that 0.32 is not well-supported is not yet published! So I simply suggest this edit; it's up to you to leave it or keep it.

Another point that would increase clarity for a general reader would be to use the terms rather than the parameter names (or, in addition to). For example, on line 210, the sentence could refer to temperature dependence and thermal niche breadth rather than just E and W0.

--------------

Reviewer #2:

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

In this article, the authors undertake several analyses to determine whether features of thermal performance curves (TPCs) in phytoplankton and prokaryotes have/can evolve or are constrained by underlying thermodynamic constraints of biochemistry. Altogether, the authors conclude that, contrary to the hypothesis that there is universal temperature dependence (UTD), that activation energy (E, describing the increase in the left side of the TPC) has diverged through time and that different clades explore considerable variation in E within clades. They also show that much of this variation has occurred in more recent parts of their timelines. They make similar conclusions that TPC breadth (Wop) displays these same patterns, which are also consistent with the evolution of respiration and photosynthesis in plants, which can be considered underlying drivers of growth. Finally, they find that E declines with latitude for all types of TPCs in their analysis, but not so for Wop.

OK, there is a lot going on in this article, so sorry it took a while. Overall, I think it's a great contribution. It provides fundamental insights into the evolution of thermal tolerance, not just retrospectively through the lens of phylogeny but also that further evolution of thermal sensitivity is possible. It's a novel take on a core problem. This paper should be the final death knell of the UTD as we currently think of it, not that there couldn't be other ways that thermodynamics constrain physiology. The paper is thorough, well-written, and the figures are engaging and informative.

I do have some questions and comments, a few little things to annoy the authors with, but hopefully they will be useful.

First, since these are most if not all (presumably) lab cultures, to what extent is adaptation/acclimation to a typical lab temperature influencing the perceived pattern of evolution? If over many generations the cultures have shifted their TPCs to adapt to ~22C, as is clearly possible given (Padfield et al. 2016), then would the results reflect some convergence around the lab environment rather than the phylogenetic pathways that these organisms travel in the wild? I don't believe there is anything to be done about this possibility in the analysis, but it could be influencing the results.

Abstract. Where it says (there are no line numbers here) 'typically-experienced thermal range', I would change this, because there is no clear reason why the rising portion of a curve is where most of the experienced temperatures are. Certainly many species can experience typical temperatures around the peak. Also see line 22 in the introduction. I would also say your results are 'consistent with' rather than 'driven by two opposing…' factors. Only if you could isolate the effect of temperature fluctuations explicitly (e.g., with experiments) could you go as far as 'driven'. Also here, evolution can occur at ecological time scales for anything; this is not a unique feature of microbes.

Line 17 - I don't think this is well-established, although I agree some folks would like it to be. See (Darveau et al. 2002).

Line 45 - I think you are overstating the restrictiveness of the UTD. Even Gillooly's original paper (Gillooly et al. 2001) stated a range of eV as "0.2 and 1.2 eV".

Line 46 - I don't think we really expect a Gaussian distribution for a measurement that cannot be negative, nor is 'very substantial' measurement error needed to generate that right skew. From my experience, even bootstrapped estimates of positive parameters are generally right-skewed without bunching up at 0. I don't think this takes away from the value of your work, but the argument seems to be a bit of a straw man. See also line 354. Maybe it's a matter of degrees.

Line 63 - Seems like you need some references here. How about (Alexander Jr and McMahon 2004).

Line 72 - Why does this require 'frequent convergence'? Seems like quickly adapting would be sufficient to break down the ability to detect heritability.

Line 141 - If you are doing this analysis phylogenetically, why are you separating the datasets at all? Can't you just run it all together?

Lines 165-168 - As I understand this here, as well as described by (Harmon et al. 2003), relative disparity varies from 0 to 1, as the subclade disparity is divided by the across-clade disparity. As described by Harmon et al, "Values near 0 imply that subclades contain relatively little of the variation present within the taxon as a whole and that, consequently, most variation is partitioned as among-subclade differences; conversely, values near 1 imply that subclades contain a substantial proportion of the total variation and thus are likely to overlap extensively". This seems to be in line with your description, but my question is that Figure 3 shows a lot of disparity above 1. Can you clarify why this can happen or confirm that there is not an error, or, clarify what I might be missing?

Line 172 - 'Exact' might be a stretch, given that the phylogenies and nodes themselves are hypotheses and estimates.

Line 176 - A bit counterintuitive here. Why is adaptive radiation characterized by decelerating evolution? Seems to me like examples of adaptive radiation are generally characterized by fast evolution accompanied by key innovations, which may or may not be part of the disparity. Why not just delete phrase after 'space'?

Line 211 - What does subtree refer to here?

Line 214 - I think optimum is the wrong word for the UTD of 0.65. As indicated in the same sentence, the value is thermodynamically-constrained, which is a biochemical/physical property and may or may not be an optimum. Likewise in line 218 and later, I would suggest a 'global mean set by thermodynamic constraints' not a 'global optimum'.

Line 261 - Might want to specify that it's a thorough analysis in prokaryotes and phytoplankton, since otherwise it suggests you have all kinds of species.

Line 267 - 'are' not 'is' here, I see two subjects.

Line 268 - Suggest here adding 'rate-limiting enzyme … of respiration'.

Line 285-286 - Not sure this argument makes sense. E doesn't really tell you breadth on its own. A low E doesn't automatically indicate narrow or wide, because the breadth is determined by the distances between upper and lower critical values and Tpk that could be connected by steep or shallow rises. I would also suggest that Wop is an arbitrary indicator of width making no assumptions about realized environmental temperatures. You might be best off getting rid of this statement, or at least getting clearer on the message.

Line 310 - Tpk in respiration, the purported underlying driver of rmax, is negatively associated with latitude (DeLong et al. 2018), which I think complicates this explanation a bit, since latitude is invoked as a proxy for fluctuations.

Line 320 - See (Uiterwaal et al. n.d.) for a recent example of a link between TPC evolution and biotic interactions as well as (Luhring and DeLong 2016) for direct evidence that biotic interactions alter rmax TPCs.

Line 348 - Not sure I see the mechanistic explanation. Even the previous paragraph explained a range of mechanisms that could underlie the changes in TPCs. Seems to me you've demonstrated patterns of evolution that clearly contradict the UTD but haven't explained the mechanisms.

Line 358-362 - I don't get this. I don't think the disparity analysis is exactly evidence of convergence, which would require demonstrating that different lineages evolved the same trait in the same conditions. It says that the clades strongly overlap in state space, which may or may not contain cases of convergence.

Signed, John P. DeLong

References

Alexander Jr, J. E., and R. F. McMahon. 2004. Respiratory response to temperature and hypoxia in the zebra mussel Dreissena polymorpha. Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology 137:425-434.

Darveau, C.-A., R. K. Suarez, R. D. Andrews, and P. W. Hochachka. 2002. Allometric cascade as a unifying principle of body mass effects on metabolism. Nature 417:166-170.

DeLong, J. P., G. Bachman, J. P. Gibert, T. M. Luhring, K. L. Montooth, A. Neyer, and B. Reed. 2018. Habitat, latitude and body mass influence the temperature dependence of metabolic rate. Biology Letters 14:20180442.

Gillooly, J. F., J. H. Brown, G. B. West, V. M. Savage, and E. L. Charnov. 2001. Effects of size and temperature on metabolic rate. Science 293:2248-2251.

Harmon, L. J., J. A. Schulte, A. Larson, and J. B. Losos. 2003. Tempo and Mode of Evolutionary Radiation in Iguanian Lizards. Science 301:961-964.

Luhring, T. M., and J. P. DeLong. 2016. Predation changes the shape of thermal performance curves for population growth rate. Current Zoology 62:501-505.

Padfield, D., G. Yvon-Durocher, A. Buckling, S. Jennings, and G. Yvon-Durocher. 2016. Rapid evolution of metabolic traits explains thermal adaptation in phytoplankton. Ecology Letters 19:133-142.

Uiterwaal, S. F., I. T. Lagerstrom, T. M. Luhring, M. E. Salsbery, and J. P. DeLong. (n.d.). Trade-offs between morphology and thermal niches mediate adaptation in response to competing selective pressures. Ecology and Evolution n/a.

--------------

Reviewer #3:

This study seeks to understand how thermal sensitivity of three key traits among prokaryotes, phytoplankton, and plants evolve along the branches of their respective phylogenetic trees. This is potentially important for knowing how species will respond to human-mediated climate change. The manuscript reports that thermal sensitivity does not evolve in a gradual manner (i.e. via random walk) but can differ considerably even between closely related species, and that there is a weak negative association between thermal sensitivity and latitude.

The paper is very nicely written - its introduction was a pleasure to read. However, in my opinion, the analyses and methods leave a great deal to be desired. The analyses are poorly described and highly unconventional (which is not a bad thing but in this case seems totally unnecessary and confused). It is not clear at all how the analyses relate to the hypotheses under investigation nor is it clear how various comparative tests should be interpreted and reconciled with each other.

I think there are some basic problems with the way the authors describe their null model - they often talk about traits evolving randomly; however, it is more accurate to say by random walk or Brownian motion. And even if a trait is evolving by random walk along the branches it does not rule out that the trait has evolved under selection - though this is a relatively minor problem. What is more serious is that the set of three analyses basically describe the same phenomenological pattern in the data but have distinct interpretations - thus the reader is left wondering how to interpret the results.

The first MCMCGLMM analyses takes all traits together and estimates heritability (h2). It is not clear how the authors constructed these multi-response models to achieve this and the performance of such models has never been assessed. Theoretically they should give identical results to more conventional tests such as Pagel's Lambda - but the prior on the random effects in the MCMCGLMM analyses can be influential. I think it would be better for the authors to conduct a more reliable test using Lambda (with or without estimating covariances among traits - depending on significance given the tree), or at least such a test should be compared with the MCMCGLMM result. In any case if we take the analyses at face value then it shows that there is some phylogenetic signal in the data, but it is not perfect - i.e. h2 is lower than one. This could be an honest indication that there is lack of signal in the trait or it could be noise in the data. The authors attempt to show that the signal is not owing to noise, but it is not clearly explained and was not convincing to me in its current state.

They then go on to try to determine if there is rate variation along and among the branches of the trees. To do this they use the 'Free model' and the 'Stable model'. These are strange choices as they don't identify significant shifts in the rate of evolution (whereas several other methods do, see below). The Free model is heavily overparameterized and at best should only be used to provide a rough visualization of rate variation (see Mooers et al, 1999 and Thomas and Freckleton, 2010 - cited in the manuscript). The Stable model was presented to estimate ancestral states rather than rates (although these are two sides of the same coin). Neither method formally identifies areas of the tree (or branches) which are evolving at a significantly different rate, yet this what the authors need to support their conclusion. Methods such as those described in Rabosky (2014, PLoS ONE 9(2): e89543), Landis et al (Systematic Biology 62(2): 193-204), Venditti et al (2011, Nature 479(7373): 393-396), Eastman et al (2011, Evolution 65(12): 3578-3589) can be used to achieve this - this list is not exhaustive . Looking at the results of the 'Free model' and the 'Stable model', most of the rate variation is at the tips of the tree which is what one would expect give the lack of phylogenetic signal in their previous analyses. It does not in itself support the idea that there are intense episodes of adaptation with closely related species having very distinct rates, which is their main conclusion.

The authors' third analyses seem superfluous. Here they estimate a model with an evolutionary optimum - this will again result in a pattern that one would expect if there was a lack of phylogenetic signal in the data. It is not clear what these analyses add to the paper (other than attractive figures) or how they should be interpreted.

There is no attempt to discriminate between the results of these first three sets of analyses which have very different interpretations despite similarities in patterns. I think the authors need to take a step back and think about how the statistical test they are carrying out relates to the hypotheses. I believe their question can be addressed by a single set of analyses to identify significant shifts in the rate of evolution on the tree using the methods I mentioned previously.

The last set of analyses regarding the latitude is also very confusing and unnecessarily complex. There is no need to estimate models with and without the variance-covariance matrix as a random effect. If there is no signal in the data, h2 will be very low - thus the variance explained by the random effect will be negligible and will reduce to be equivalent to the non-phylogenetic test.

Owing to the reasons outlined above, while I find the subject matter interesting and potentially suitable for PLoS Biology, I can't support its publication in the current form. In addition, based on the analyses provided I cannot predict whether the result will be of interest to the readership of PLoS Biology.

3 Apr 2020

Submitted filename: Kontopoulos_etal_Rev_Responses.pdf

Click here for additional data file.

11 May 2020

Dear Dr Kontopoulos,

Thank you very much for submitting a revised version of your manuscript "Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate" for consideration as a Research Article at PLOS Biology. This revised version of your manuscript has been evaluated by the PLOS Biology editors, the Academic Editor and two of the original reviewers.

In light of the reviews (below), we offer you the opportunity to address the remaining points from the reviewers in a revised version that we anticipate should not take you very long. We will then assess your revised manuscript and your response to the reviewers' comments and we may consult the reviewers again. I should warn you, however, that we will only consult the reviewers one more time, and if they remain unsatisfied then we will not consider the paper further.

We expect to receive your revised manuscript within 1 month.

Please email us (plosbiology@plos.org) if you have any questions or concerns, or would like to request an extension. At this stage, your manuscript remains formally under active consideration at our journal; please notify us by email if you do not intend to submit a revision so that we may end consideration of the manuscript at PLOS Biology.

**IMPORTANT - SUBMITTING YOUR REVISION**

Your revisions should address the specific points made by each reviewer. Please submit the following files along with your revised manuscript:

1. A 'Response to Reviewers' file - this should detail your responses to the editorial requests, present a point-by-point response to all of the reviewers' comments, and indicate the changes made to the manuscript.

*NOTE: In your point by point response to the reviewers, please provide the full context of each review. Do not selectively quote paragraphs or sentences to reply to. The entire set of reviewer comments should be present in full and each specific point should be responded to individually.

You should also cite any additional relevant literature that has been published since the original submission and mention any additional citations in your response.

2. In addition to a clean copy of the manuscript, please also upload a 'track-changes' version of your manuscript that specifies the edits made. This should be uploaded as a "Related" file type.

*Resubmission Checklist*

When you are ready to resubmit your revised manuscript, please refer to this resubmission checklist: https://plos.io/Biology_Checklist

To submit a revised version of your manuscript, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' where you will find your submission record.

Please make sure to read the following important policies and guidelines while preparing your revision:

*Published Peer Review*

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*PLOS Data Policy*

Please note that as a condition of publication PLOS' data policy (http://journals.plos.org/plosbiology/s/data-availability) requires that you make available all data used to draw the conclusions arrived at in your manuscript. If you have not already done so, you must include any data used in your manuscript either in appropriate repositories, within the body of the manuscript, or as supporting information (N.B. this includes any numerical values that were used to generate graphs, histograms etc.). For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5

*Blot and Gel Data Policy*

We require the original, uncropped and minimally adjusted images supporting all blot and gel results reported in an article's figures or Supporting Information files. We will require these files before a manuscript can be accepted so please prepare them now, if you have not already uploaded them. Please carefully read our guidelines for how to prepare and upload this data: https://journals.plos.org/plosbiology/s/figures#loc-blot-and-gel-reporting-requirements

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosbiology/s/submission-guidelines#loc-materials-and-methods

Thank you again for your submission to our journal. We hope that our editorial process has been constructive thus far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Roli Roberts

Roland G Roberts, PhD,

Senior Editor

PLOS Biology

*****************************************************

REVIEWERS' COMMENTS:

Reviewer #2:

[and see attached figure]

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

I previously reviewed this article. Overall, I think the paper is improved and the responses were reasonable. I just have a few areas that I still find challenging.

1) The disparity levels going above 1 still seems confusing. The authors write "The disparity-through-time method compares i) the mean disparity within subclades with ii) the total disparity across the phylogeny. Therefore, it is mathematically possible to get values much greater than 1 when the trait variance within subclades (at a given time point) is much greater than the trait variance across the tree." What I don't understand is how there can be more variance within than across subclades, when the variation across subclades would seem to me to be inclusive of all the variation occurring within subclades. I suspect there is some crucial detail missing.

2) The use of Wop and the link between E and breadth is still confusing. There are two issues.

2A) First, choosing Wop as you have done does not require an assumption. It is just a choice. I think what you might be trying to say is that you chose this specific metric (instead of, say, the difference between the upper and lower critical temperatures) because you feel this metric best captures the curve's breadth in the typical range of experienced temperatures.

2B) Second, the explanation that E reflects breadth still doesn't make sense to me. The authors write "an increase in Wop (a larger operational niche width) will necessarily be associated with a shallower rising slope (a low E value)." Yet it is easy to show that this is not true. Consider the attached cartoon. Two curves, same Tpk. The lower curve has a lower E and a smaller Wop - opposite of the 'necessary' pattern. So I find this confusing, and it influences interpretation of the results as well. For example, the explanation for Figure S4 is that it has to be that way. Does it? I think, perhaps, the authors are making additional unstated assumptions about correlations among parameters, such that the proposed necessary correlation will arise under particular scenarios, such as, for example, Tpk and E are negatively correlated.

3) It is still not very clear why the disparity analysis indicates convergent evolution. The authors clarify that there are different definitions of convergence that focus on pattern or process, and that they 'use the term "convergent evolution" to refer to remote clades that increasingly overlap in the parameter space due to independent evolution of similar thermal sensitivity values.' Note the point here is 'similar thermal sensitivity values'. To me the disparity results seem to say that the subclades have evolved similar 'ranges' - that is, a lot of overlap in extent - rather than showing convergence toward similar values within that extent. The pattern appears to be (and this is a major point of the paper) that lineages are exploring a lot of thermal parameter space. Overlapping a lot because multiple lineages have evolved a wide range of parameter values does not seem to me to be consistent with the idea of convergent evolution.

A few more minor things:

Figure 1 legend - "E and ED (eV) control how smoothly". Isn't it how steeply, since a low E does not cause a 'rough' rise?

L123 - Evolution would not favor species but individuals with lower thermal sensitivities.

L214 - The 'true' disparity is not really knowable. The 'estimated' or 'observed' disparity is - 'observed' is the word you used in the F3 legend. I would use that.

Reviewer #3:

I will speak to the author's comments using the numbers they include in their rebuttal:

30) I think the removal of the notion of gradual evolution has improved the manuscript.

31) The authors are misguided in their assessment: the MCMCglmm analyses are not far superior. h2 and Lambda are mathematically identical (if for lambda you take a mean of the within-species values and estimate co-variances among traits) - this is clear in the original paper describing heritability in the context of phylogenetic analysis, which the authors cite for a description of the model! My specific worry is that in the MCMCglmm implementation, h2 can be very sensitive to the prior. (Incidentally, what prior was on the random effect here, and what does the distribution look like if you turn the likelihood function off and just target the prior?) This sensitivity to prior choice is not the case for Lambda, owing to the fact that it is implemented within the more conventional GLS framework. I therefore still maintain that given this, it is important to check the test using the more well-established method. Estimating co-variances and Lambda can be done in many GLS packages including BayesTraits - Pagel's own program. Incidentally, phylogenetic imputation can also be done in a GLS framework again in number of packages including BayesTraits.

32) Given that the Levy process analyses have been conducted and are by far better justified I believe they should be reported in the main text - the others can be removed altogether as they are no longer needed.

33) I still can't see the value of these analyses if authors believe the previous analyses - they all show the same thing - but the interpretations are not consistent. I believe they must be removed - how can these results and the Levy process results be reconciled?

34) The reasons I state above mean we cannot believe the h2 until it has been verified by lambda so this point still stands.

Submitted filename: Presentation1.jpg

Click here for additional data file.

9 Jun 2020

Submitted filename: responses_to_the_Reviewers.pdf

Click here for additional data file.

18 Jul 2020

Dear Dr Kontopoulos,

Thank you very much for submitting a revised version of your manuscript "Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate" for consideration as a Research Article at PLOS Biology. This revised version of your manuscript has been evaluated by the PLOS Biology editors, the Academic Editor and one of the original reviewers.

In light of the reviews (below), we are pleased to offer you the opportunity to address the remaining points from reviewer #3 in a revised version that we anticipate should not take you very long.

We note that we had previously said that we would not consider the manuscript further if you failed to satisfy the reviewers; however, because this reviewer remains overall very positive about your study, we are prepared to give you one final chance. Further consideration is absolutely predicated on you addressing this final issue. We will then assess your revised manuscript and your response to the reviewer's comments and we may consult this reviewer again.

We expect to receive your revised manuscript within 1 month.

Please email us (plosbiology@plos.org) if you have any questions or concerns, or would like to request an extension. At this stage, your manuscript remains formally under active consideration at our journal; please notify us by email if you do not intend to submit a revision so that we may end consideration of the manuscript at PLOS Biology.

**IMPORTANT - SUBMITTING YOUR REVISION**

Your revisions should address the specific points made by each reviewer. Please submit the following files along with your revised manuscript:

1. A 'Response to Reviewers' file - this should detail your responses to the editorial requests, present a point-by-point response to all of the reviewers' comments, and indicate the changes made to the manuscript.

*NOTE: In your point by point response to the reviewers, please provide the full context of each review. Do not selectively quote paragraphs or sentences to reply to. The entire set of reviewer comments should be present in full and each specific point should be responded to individually.

You should also cite any additional relevant literature that has been published since the original submission and mention any additional citations in your response.

2. In addition to a clean copy of the manuscript, please also upload a 'track-changes' version of your manuscript that specifies the edits made. This should be uploaded as a "Related" file type.

*Resubmission Checklist*

When you are ready to resubmit your revised manuscript, please refer to this resubmission checklist: https://plos.io/Biology_Checklist

To submit a revised version of your manuscript, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' where you will find your submission record.

Please make sure to read the following important policies and guidelines while preparing your revision:

*Published Peer Review*

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*PLOS Data Policy*

Please note that as a condition of publication PLOS' data policy (http://journals.plos.org/plosbiology/s/data-availability) requires that you make available all data used to draw the conclusions arrived at in your manuscript. If you have not already done so, you must include any data used in your manuscript either in appropriate repositories, within the body of the manuscript, or as supporting information (N.B. this includes any numerical values that were used to generate graphs, histograms etc.). For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5

*Blot and Gel Data Policy*

We require the original, uncropped and minimally adjusted images supporting all blot and gel results reported in an article's figures or Supporting Information files. We will require these files before a manuscript can be accepted so please prepare them now, if you have not already uploaded them. Please carefully read our guidelines for how to prepare and upload this data: https://journals.plos.org/plosbiology/s/figures#loc-blot-and-gel-reporting-requirements

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosbiology/s/submission-guidelines#loc-materials-and-methods

Thank you again for your submission to our journal. We hope that our editorial process has been constructive thus far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Roli Roberts

Roland G Roberts, PhD,

Senior Editor

PLOS Biology

*****************************************************

REVIEWER'S COMMENTS:

Reviewer #3:

I find the difference between Lambda and H2 quite alarming - and I do not really understand why the authors would not use the MCMC option in a Bayesian context for comparison which allows multiple data points per taxa - especially as MCMCglmm also provides a posterior distribution of h2. It would be ideal to see the two distributions together. Such differences in the level of phylogenetic signal have the potential to change the results altogether - I am not sure if that is the case here as I don't have enough information. I think this is a critical point as phylogenetic signal is fundamental to this study. I feel like I am sounding like a broken record here, but the analyses I am suggesting a very easy to do and I can't understand why the authors don't want to perform/show them.

I very much like this paper and hope to see it published on PLoS Biology - but I think this final piece of the puzzle is missing.

28 Jul 2020

Submitted filename: response_to_the_reviewer.pdf

Click here for additional data file.

12 Aug 2020

Dear Dr Kontopoulos,

Thank you for submitting your revised Research Article entitled "Adaptive evolution shapes the present-day distribution of the thermal sensitivity of population growth rate" for publication in PLOS Biology. I have now obtained further advice from reviewer #3 and have discussed their comments with the Academic Editor.

You'll see that reviewer #3 is now largely satisfied, and the Academic Editor agrees that this reviewer's final point can be addressed textually. Based on the review, we will probably accept this manuscript for publication, assuming that you will modify the manuscript to address the remaining points raised by reviewer #3. Please also make sure to address the data and other policy-related requests noted at the end of this email.

We expect to receive your revised manuscript within two weeks. Your revisions should address the specific points made by each reviewer. In addition to the remaining revisions and before we will be able to formally accept your manuscript and consider it "in press", we also need to ensure that your article conforms to our guidelines. A member of our team will be in touch shortly with a set of requests. As we can't proceed until these requirements are met, your swift response will help prevent delays to publication.

*Copyediting*

Upon acceptance of your article, your final files will be copyedited and typeset into the final PDF. While you will have an opportunity to review these files as proofs, PLOS will only permit corrections to spelling or significant scientific errors. Therefore, please take this final revision time to assess and make any remaining major changes to your manuscript.

NOTE: If Supporting Information files are included with your article, note that these are not copyedited and will be published as they are submitted. Please ensure that these files are legible and of high quality (at least 300 dpi) in an easily accessible file format. For this reason, please be aware that any references listed in an SI file will not be indexed. For more information, see our Supporting Information guidelines:

https://journals.plos.org/plosbiology/s/supporting-information

*Published Peer Review History*

Please note that you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out. Please see here for more details:

https://blogs.plos.org/plos/2019/05/plos-journals-now-open-for-published-peer-review/

*Early Version*

Please note that an uncorrected proof of your manuscript will be published online ahead of the final version, unless you opted out when submitting your manuscript. If, for any reason, you do not want an earlier version of your manuscript published online, uncheck the box. Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us as soon as possible if you or your institution is planning to press release the article.

*Protocols deposition*

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosbiology/s/submission-guidelines#loc-materials-and-methods

*Submitting Your Revision*

To submit your revision, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' to find your submission record. Your revised submission must include a cover letter, a Response to Reviewers file that provides a detailed response to the reviewers' comments (if applicable), and a track-changes file indicating any changes that you have made to the manuscript.

Please do not hesitate to contact me should you have any questions.

Sincerely,

Roli Roberts

Roland G Roberts, PhD,

Senior Editor,

rroberts@plos.org,

PLOS Biology

------------------------------------------------------------------------

DATA POLICY:

You may be aware of the PLOS Data Policy, which requires that all data be made available without restriction: http://journals.plos.org/plosbiology/s/data-availability. For more information, please also see this editorial: http://dx.doi.org/10.1371/journal.pbio.1001797

We note that you are intending to deposit your phylogeny and parameter estimates in Figshare; please could you do so, include the URL in the manuscript, and provide us with access to it? Please could you also deposit (e.g. in Github) any code used for the study?

In addition, we ask that all individual quantitative observations that underlie the data summarized in the figures and results of your paper be made available in one of the following forms:

1) Supplementary files (e.g., excel). Please ensure that all data files are uploaded as 'Supporting Information' and are invariably referred to (in the manuscript, figure legends, and the Description field when uploading your files) using the following format verbatim: S1 Data, S2 Data, etc. Multiple panels of a single or even several figures can be included as multiple sheets in one excel file that is saved using exactly the following convention: S1_Data.xlsx (using an underscore).

2) Deposition in a publicly available repository. Please also provide the accession code or a reviewer link so that we may view your data before publication.

Regardless of the method selected, please ensure that you provide the individual numerical values that underlie the summary data displayed in all of the main and Supplementary figure panels as they are essential for readers to assess your analysis and to reproduce it. NOTE: the numerical data provided should include all replicates AND the way in which the plotted mean and errors were derived (it should not present only the mean/average values).

Please also ensure that figure legends in your manuscript include information on where the underlying data can be found, and ensure your supplemental data file/s has a legend.

Please ensure that your Data Statement in the submission system accurately describes where your data can be found.

------------------------------------------------------------------------

REVIEWER'S COMMENTS:

Reviewer #3:

One can see that the results of the phylogenetic signal are quite different from different methods! As the authors are happy to include this and discuss them I am happy to recommend publication – with one caveat. The authors seem to think that the difference is associated with the way the program deal with missing data. I doubt very much this is the case – one should not be able to get more information about signal from imputed/missing data points! I am sure that the difference between mcmcglmm and Baystraits is the prior. The prior is far more simple to implement in a GLS framework (BayesTraits) compared to the way it has to be implemented in mcmcglmm. I think the authors should at least posit this idea – then the reader can decide and/or test it in the future.

1 Sep 2020

Submitted filename: response_to_the_reviewer.pdf

Click here for additional data file.

14 Sep 2020

Dear Dr Kontopoulos,

On behalf of my colleagues and the Academic Editor, Simon A. Levin, I am pleased to inform you that we will be delighted to publish your Research Article in PLOS Biology.

The files will now enter our production system. You will receive a copyedited version of the manuscript, along with your figures for a final review. You will be given two business days to review and approve the copyedit. Then, within a week, you will receive a PDF proof of your typeset article. You will have two days to review the PDF and make any final corrections. If there is a chance that you'll be unavailable during the copy editing/proof review period, please provide us with contact details of one of the other authors whom you nominate to handle these stages on your behalf. This will ensure that any requested corrections reach the production department in time for publication.

Early Version

The version of your manuscript submitted at the copyedit stage will be posted online ahead of the final proof version, unless you have already opted out of the process. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

PRESS

We frequently collaborate with press offices. If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximise its impact. If the press office is planning to promote your findings, we would be grateful if they could coordinate with biologypress@plos.org. If you have not yet opted out of the early version process, we ask that you notify us immediately of any press plans so that we may do so on your behalf.

We also ask that you take this opportunity to read our Embargo Policy regarding the discussion, promotion and media coverage of work that is yet to be published by PLOS. As your manuscript is not yet published, it is bound by the conditions of our Embargo Policy. Please be aware that this policy is in place both to ensure that any press coverage of your article is fully substantiated and to provide a direct link between such coverage and the published work. For full details of our Embargo Policy, please visit http://www.plos.org/about/media-inquiries/embargo-policy/.

Thank you again for submitting your manuscript to PLOS Biology and for your support of Open Access publishing. Please do not hesitate to contact me if I can provide any assistance during the production process.

Kind regards,

Vita Usova,

Publishing Editor

PLOS Biology

on behalf of

Roland Roberts,

Senior Editor

PLOS Biology