Author response: Latent functional diversity may accelerate microbial community responses to temperature fluctuations

Abstract

Article Figures and data Abstract Editor's evaluation eLife digest Introduction Materials and methods Results Discussion Data availability References Decision letter Author response Article and author information Metrics Abstract How complex microbial communities respond to climatic fluctuations remains an open question. Due to their relatively short generation times and high functional diversity, microbial populations harbor great potential to respond as a community through a combination of strain-level phenotypic plasticity, adaptation, and species sorting. However, the relative importance of these mechanisms remains unclear. We conducted a laboratory experiment to investigate the degree to which bacterial communities can respond to changes in environmental temperature through a combination of phenotypic plasticity and species sorting alone. We grew replicate soil communities from a single location at six temperatures between 4°C and 50°C. We found that phylogenetically and functionally distinct communities emerge at each of these temperatures, with K-strategist taxa favored under cooler conditions and r-strategist taxa under warmer conditions. We show that this dynamic emergence of distinct communities across a wide range of temperatures (in essence, community-level adaptation) is driven by the resuscitation of latent functional diversity: the parent community harbors multiple strains pre-adapted to different temperatures that are able to ‘switch on’ at their preferred temperature without immigration or adaptation. Our findings suggest that microbial community function in nature is likely to respond rapidly to climatic temperature fluctuations through shifts in species composition by resuscitation of latent functional diversity. Editor's evaluation This important study tests potential mechanisms for microbial community adaptation to temperature. Using elegant experiments, the authors convincingly show that selection on standing variation in the community, that is, species sorting, drives the community response to temperature. This article will be of interest to ecologists and microbiologists studying the impacts of global change on community and ecosystem processes. https://doi.org/10.7554/eLife.80867.sa0 Decision letter Reviews on Sciety eLife's review process eLife digest Most ecosystems on Earth rely on dynamic communities of microorganisms which help to cycle nutrients in the environment. There is increasing concern that climate change may have a profound impact on these complex networks formed of large numbers of microbial species linked by intricate biochemical relationships. Any species within a microbial community can acclimate to new temperatures by quickly tweaking their biological processes, for example by activating genes that are more suited to warmer conditions. Over time, a species may acclimate or adapt to new conditions. However, the community as a whole can also respond to these changes, and often much faster, by simply altering the abundance or presence of its members through a process known as species sorting. It remains unclear exactly how acclimation, adaptation and species sorting each contribute to the community’s response to a temperature shift – an increasingly common scenario under global climate change. To address this question, Smith et al. investigated how species sorting and acclimation may help whole soil bacterial communities to cope with lasting changes in temperature. To do so, soil samples from a single field site (and therefore featuring the same microbial community) were incubated for four weeks under six different temperatures. Genetic analyses revealed that, at the end of the experiments, distinct communities specific to a given temperature had emerged. They all differed in species composition and the types of biological functions they could perform. Further experiments showed that each community had been taken over by strains of bacteria which grew best at the new temperature that they had been exposed to, including extreme warming scenarios never seen in their native environment. This suggests that these organisms were already present in the original community. They had persisted even under temperatures which were not optimal for them, acting as a slumbering (‘latent’) ‘reservoir’ of traits and functional abilities that allowed species sorting to produce distinct and functionally capable communities in each novel thermal environment. This suggests that species sorting could help bacterial communities to cope with dramatic changes in their thermal environment. Smith et al.’s findings suggest that bacterial communities can cope with warming environments much better than has been previously thought. In the future, this work may help researchers to better predict how climate change could impact microbial community structure and functioning, and most crucially their contributions to the global carbon cycle. Introduction Microbes are drivers of key ecosystem processes. They are tightly linked to the wider ecosystem as pathogens, mutualists, and food sources for higher trophic levels, and also play a central role in ecosystem-level nutrient cycling, and therefore, ultimately in global biogeochemical cycles. Temperature has a pervasive influence on microbial communities because of its direct impact on microbial physiology and fitness (Oliverio et al., 2017; García et al., 2018; Smith et al., 2019). There is therefore great interest in understanding how temperature fluctuations impact microbial community dynamics and how those impacts affect the wider ecosystem (Bardgett et al., 2008). Temperature varies at practically all biologically relevant timescales, from seconds (e.g., sun/shade), through daily and seasonal fluctuations, to longer-term changes, including anthropogenic climate warming and fluctuations over geological timescales. Whole microbial communities can respond to temperature changes over time and space through phenotypic (especially, physiological) plasticity (henceforth, ‘acclimation’), as well as genetic adaptation in their component populations (Bennett et al., 1990; Kishimoto et al., 2010; Blaby et al., 2012; Kontopoulos et al., 2020a). Microbial thermal acclimation can occur relatively rapidly (timescales of minutes to days) through processes such as activation and up- or downregulation of particular genes and alteration of fatty acids used in building cell walls (Suutari and Laakso, 1994). Adaptation is a necessarily slower process (timescales of weeks or longer) occurring either through selection on standing genetic variation in the population or that arising through recombination and mutation (Bennett et al., 1990; Padfield et al., 2016; Barton et al., 2020). In addition, a third key mechanism through which microbial communities can respond to changing temperatures is species sorting (Leibold et al., 2004; Wu et al., 2018): changes in community composition through species-level selection where taxa maladapted to a new temperature are replaced by those that are pre-adapted to it. This can happen either relatively rapidly through the resuscitation or suppression of taxa that are already present (Lennon and Jones, 2011; Wisnoski and Lennon, 2021), or more slowly through immigration-extinction dynamics driven by dispersal from the regional species pool (Langenheder and Székely, 2011; Wu et al., 2018). Resuscitation may be an important mechanism driving species sorting in microbial communities in particular because many microbial taxa have the capacity to form environment-resistant spores when conditions are unfavorable, and then rapidly activate metabolic pathways and resuscitate in favorable conditions. This effectively widens their thermal niche to allow persistence in the face of temperature change (Lennon and Jones, 2011; Wisnoski and Lennon, 2021). In order for rapid resuscitation of dormant taxa to allow species sorting to drive community-level adaptation, there must be a wide source pool of species to select from. Indeed, sequencing studies have revealed the presence of thousands of distinct microbial taxa in small environmental samples, most occurring at low abundance (Lynch and Neufeld, 2015; Sogin et al., 2006; Thompson et al., 2017). There is also strong evidence that bacteria are often found well outside of their thermal niche. For example, thermophilic taxa are often found in cold ocean beds and cool soils (Marchant et al., 2008; Hubert et al., 2009; Zeigler, 2014). Thus, a significant reservoir of latent microbial functional diversity may be commonly present for species sorting to act upon (Lennon and Jones, 2011; Wisnoski and Lennon, 2021). Understanding the relative importance of acclimation, adaptation, and species sorting in the assembly and turnover (succession) of microbial communities is key to determining the rate at which they can respond to different regimes of temperature fluctuations. For example, a combination of acclimation and species sorting through resuscitation would enable rapid responses to sudden temperature changes, relative to adaptation. A number of past studies have investigated responses of microbial community composition and functioning to temperature changes, showing that composition can respond rapidly to warming (Allison and Martiny, 2008; Aydogan et al., 2018), often correlated with responses of ecosystem functioning (Karhu et al., 2014; Melillo et al., 2017; Yu et al., 2018). However, a mechanistic basis of these community-level responses remains elusive, both in terms of how individual taxa respond to changing temperatures in a community context and the relative importance of acclimation, adaptation, and species sorting. The community context of the responses of individual microbial populations is important because interactions between strains can constrain or accelerate acclimation as well as adaptive evolution (Scheuerl et al., 2020). Also, while the importance of species sorting in microbial communities per se has been studied (Van der Gucht et al., 2007; Langenheder and Székely, 2011; Székely and Langenheder, 2014), work on this issue in the context of environmental temperature is practically nonexistent. A further consideration is whether differing temperature conditions, such as the frequency and magnitude of temperature fluctuations, may influence the life history strategies of the taxa in the community (Gilchrist, 1995; Basan et al., 2020), which will in turn alter the relative importance of sorting, acclimation, and adaption. In order to identify the life history strategies of bacteria, we must quantify their phenotypic traits, such as growth rates and yield (Malik et al., 2020). Quantifying these traits can allow us to identify growth specialists (r-strategists) and carrying-capacity specialists (K-strategists) (Marshall, 1986), and thus test whether these strategies are differentially favored in different thermal environments. By identifying life history strategies, we can consider the ecosystem implications of any adaptation-, acclimation-, or sorting-driven changes in microbial communities (Malik et al., 2020). Here, we investigate whether species sorting and latent functional diversity alone can influence the response of soil bacterial communities to changes in environmental temperature. To this end, we subject replicate communities, shielded from immigration, to a wide range of temperatures in the laboratory. In order to understand the mechanistic basis of observed community-level changes, we analyze the phylogenetic structure and functional traits of the resulting component taxa. Materials and methods We performed a species-sorting experiment to investigate how microbial communities respond to shifts in temperature (Figure 1). After each community incubation at a given temperature, we estimated the thermal optimum (Topt) for every isolated strain by measuring the thermal performance curve (TPC) of its maximal growth rate across several temperatures (Figure 1D). This allowed us to determine how strain-level thermal preferences and niche widths vary with community growth (isolation) temperature, and the presence of taxa pre-adapted to the new temperature. We also performed a phylogenetic analysis of the overall assemblage to identify whether deep evolutionary differences predict which taxa (and their associated traits) are favored by sorting at different temperatures. To quantify strain-level functional traits, we measured their available cellular metabolic energy (ATP), respiration rates, and biomass yield at population steady state (carrying capacity), which allowed us to identify r- vs. K-strategists as well as trade-offs between different strategies. Figure 1 Download asset Open asset The species sorting experiment. (A) Different bacterial taxa (colored circles) sampled from the soil community. (B) Samples maintained at 4, 10, 21, 30, 40, and 50°C (only three temperatures shown for illustration), allowing species sorting for 4 weeks. (C) Soil washes from each core plated out onto agar and grown at both the sorting temperature and 22°C (standard temperature) to allow further species sorting and facilitate isolation (next step). (D) The six most abundant (morphologically different) colonies from each plate were picked, streaked, and isolated, and their physiological and life history traits measured. The curves represent each strain’s unique unimodal response of growth rate to temperature. Species-sorting experiment Soil cores were taken from a single site in Nash’s Field (Silwood Park, Berkshire, UK, the site of a long-term field experiment [Macdonald et al., 2015]) in June 2016 (Figure 1A). Six cores were taken from the top 10 cm of soil, using a 1.5-cm-diameter sterile corer. Ambient soil temperature at the time of sampling was 19.4°C. The cores were maintained at different temperatures in the laboratory (4, 10, 21, 30, 40, and 50°C) for 4 weeks to allow species sorting to occur at those temperatures (Figure 1B). The soil was rehydrated periodically with sterile, deionized water during incubation. During this period, in each microcosm (incubated soil core), we expected some taxa would go extinct if the temperature was outside their thermal niche, and that survivors would acclimate to the new local thermal conditions. We also expected that the 4-week incubation period would be sufficient time for changes to species interactions due to changes in abundance or traits, and therefore that interaction-driven sorting would occur in addition to the immediate extinctions and acclimation. Because bacteria display higher growth rates at warmer temperatures (Smith et al., 2019), the different incubation conditions could result in differential generational turnover of species across the given timescale. However, we did not supplement the soil samples with any additional nutrients and thus expect any growth of bacteria during this time to be heavily restricted due to nutrient limitation. Therefore, environmental exclusion (elimination of taxa maladapted to the temperature conditions) was expected to be the dominant process affecting the bacterial taxa during this stage of the sorting experiment, rather than changes in abundances due to population growth. We then isolated bacterial strains by washing the soil with PBS, plating the soil wash onto R2 agar, and incubating the plates at both their 4-week incubation temperature treatments (‘sorting temperature’) and at 22°C (‘standard temperature’). The sorting temperature allowed us to determine whether strains in each community tended to have thermal optima-matching experimental temperatures, while the standard temperature allowed us to determine whether a 4-week incubation resulted in a loss of taxa that were poorly adapted to 22°C. Appearance of strains with thermal optima matching the standard temperature would indicate incomplete species sorting because the 4-week treatment at temperatures higher or lower than 22°C had not eliminated (or at least suppressed) them. The plates were incubated until bacterial colonies formed, of which we isolated a single colony from each of the six most abundant morphologically distinct colony types on each plate (Figure 1C). Additional species sorting likely occurred during this plating-based isolation because strains with the highest growth rates at each temperature would be the first to form visible colonies and be selected. The time frame for colony appearance on the agar plates differed between temperature treatments, ranging from (∼10 days at 4°C to ∼1.5 days at 50°C). Morphologically distinct colonies were isolated from each of the six sorting-temperature and six standard-temperature plates on R2 agar by streak-plating, before being frozen as glycerol stocks (Figure 1), which were later revived for trait measurements (see below). In total, 74 strains were isolated in this way. Taxonomic identification Request a detailed protocol 16S rDNA sequences were used to identify the isolates. Raw sequences were first trimmed using Geneious 10.2.2 (https://www.geneious.com), and BLAST searches were then used to assign taxonomy to each trimmed sequence at the genus level. GenBank accession numbers of sequences are provided in Table 2. Quantifying physiological and life history traits Growth, respiration, and ATP content Request a detailed protocol We measured growth rate and respiration rate simultaneously across a range of temperatures for each isolate to construct its acute TPCs for these two traits. We henceforth denote the maximum growth rate across the temperature range by μmax, and the temperature at which this growth rate maximum occurs as Topt (optimal growth temperature or thermal optimum). The ATP content of the entire cell culture was also measured at the start and end of the growth assay. Strains were revived from glycerol stocks into fresh LB broth and incubated to carrying capacity at the temperature of the subsequent experiment. This growth to carrying capacity was an acclimation period, which typically took between 72 hr (warmest temperatures) to 500 hr (coldest temperature). Biomass abundance was determined by OD600 – optical density measurements at 600 nm wavelength. Prior to each growth-respiration assay, the strains were diluted 1:100 in LB, pushing them into a short lag phase before exponential growth started again (also tracked by OD600 measurements). The exponentially growing cultures were subsequently centrifuged at 8000 rpm for 5 min to pellet the cells, which were then resuspended in fresh LB to obtain 400 µl culture at a final OD600 of ∼0.2–0.3. This yielded cells primed for immediate exponential growth without a lag phase. These cultures were serially diluted in LB (50% dilutions) three times, producing a range of starting densities of growing cells (four biological replicates per strain/temperature combination). 100 µl subcultures of each replicate population were taken and OD600 was tracked in a Synergy HT microplate reader (BioTek Instruments, USA) to ensure that cells were indeed in exponential growth. Initial ATP measurements were made using the BacTiter-Glo assay (see below for details) and cell counts were taken using a BD Accuri C6 flow cytometer (BD Biosciences, USA). Cells were then incubated with a MicroResp plate to capture cumulative respiration (see below for details of the MicroResp system) at the experimental temperature and allowed to continue growing for a short period of time (typically 3–4 hr). After growth, the MicroResp plate was read, and final cell count and ATP measurements taken. We estimated average cell volumes and calculated the cellular carbon per cell from the flow cytometry cell diameter measurements using the relationship from Romanova and Sazhin, 2010: fgC cell-1=133.754⁢V0.438. Multiplying this by the cell counts gives an estimate of the carbon biomass of the culture at the starting and ending points. The difference between the initial biomass and biomass at the end of the experiment gives the total carbon sequestered through growth. Given an initial biomass (C0) that grows over time (t) to reach a final biomass (Ctot), assuming the population is in exponential growth, the mass-specific growth rate (µ) is given by μ=log⁡(CtotC0)t. Respiration rates of cultures were measured during growth using the MicroResp system (Campbell et al., 2003). This is a colorimetric assay initially developed to measure CO2 production from soil samples, which has since been used to measure respiration of bacterial cultures (Lawrence et al., 2012; Foster and Bell, 2012; Rivett et al., 2017). We calculate the biomass-specific respiration rate using an equation that accounts for changes in biomass of the growing cultures over time (Smith et al., 2021): R=μ⁢RtotC0⁢eμ⁢t-C0. Here, Rtot is the total mass of carbon produced according to the MicroResp measurements, C0 is the initial population biomass, µ is the previously calculated growth rate, and t is the experiment duration. ATP content of the cultures was measured using the Promega BacTiter-Glo reagent, which produces luminescence in the presence of ATP, proportional to the concentration of ATP. 50 µl of culture (diluted 1:100) was incubated with 25 µl reagent. Luminescence was measured over a 6 min period to allow the reaction to develop completely, and measurements of luminescence recorded at the 0, 2, 4, and 6 min timepoints. The highest relative light unit (RLU) measurement for each culture was taken and used to calculate the quantity of ATP, using log⁡(nM ATP)=1.21⋅log⁡(RLU)-4.69, derived from a calibration curve. This was then normalized by the flow cytometry measurements to calculate the value of ATP/biomass. Thermal performance curves Request a detailed protocol To quantify TPCs of individual isolates, we fitted the Sharpe–Schoolfield model with the temperature of peak performance (Tpk) as an explicit parameter (Schoolfield et al., 1981; Kontopoulos et al., 2020b) to the experimentally derived temperature-dependent growth rate and respiration rates of each isolate: (1) B⁢(T)=B0⁢e-Ek⋅(1T-1Tref)1+EED-E⁢eEDk⁢(1Tpk-1T). Here, T is the temperature in Kelvin (K), B is the biological rate (in this case, either growth rate, µ, or respiration rate, R), B0 is the temperature-independent metabolic rate constant approximated at some (low) reference temperature Tref, E is the activation energy in electron volts (eV) (a measure of ‘thermal sensitivity’), k is the Boltzmann constant (8.617×10-5 eV K-1), Tpk is the temperature where the rate peaks, and ED is the deactivation energy, which determines the rate of decline in the biological rate beyond Tpk. We then calculated the peak performance (i.e., Rmax or μmax) by solving Equation 1 for T=Tpk. This model was fitted to each dataset using a standard nonlinear least-squares procedure (Smith et al., 2021). The Tpk for growth rate was considered the optimum growth temperature (i.e., Topt) for each isolate. Then, the operational niche width was calculated as the difference between Topt and the temperature below this value where μmax (B⁢(T) in Equation 1) reached 50% of its maximum (i.e., μmax at Topt). This, a measure of an organism’s thermal niche width relevant to typically experienced temperatures (Pawar et al., 2016; Kontopoulos et al., 2020a), was used as a quantification of the degree to which an isolate is a thermal generalist or specialist. In most cases, Topt was derived directly from the Sharpe–Schoolfield flow cytometry growth rate fits. Four strains of Streptomyces were unsuitable for standard flow cytometry methods due to their formation of mycelial pellets (van Veluw et al., 2012). For these strains, growth rates derived from optical density measurements were used to estimate Topt instead. Trade-offs between traits Request a detailed protocol To understand the trade-offs and collinearities between different life history and physiological traits, we performed a principal components analysis (PCA), with optimum growth temperature (Topt), niche width, peak growth rate (μmax), peak respiration rate (Rmax), mean cellular ATP content (log-transformed), and carrying capacity (OD600) as input variables (scaled to have mean = 0 and SD = 1). All rate calculations, model fitting, and analyses were performed in R. Comparison to alternative datasets Request a detailed protocol We additionally investigated phylum-level life history strategy differences in two previously collated meta-analysis datasets as a comparison to our findings. DeLong et al., 2010 compiled data on both active (growth phase) and passive (stationary phase) metabolic rates, as well as growth rates, across a range of bacteria (mainly from Makarieva et al., 2005), which were corrected to 20°C using an activation energy of 0.61 eV. We also investigated differences in the growth rates of bacteria compiled in Smith et al., 2019, which we temperature-corrected to 20°C here for comparison to the DeLong et al., 2010 dataset, based on each strain’s individual TPC parameters. Phylogenetic trait mapping Request a detailed protocol We used 16S sequences to build a phylogeny in order to investigate the evolution of thermal performance across the isolated bacterial taxa. Sequences were aligned to the SILVA 16S reference database using the SILVA Incremental Aligner (SINA) (Pruesse et al., 2012). From this alignment, 100 trees were inferred in RAxML (v8.1.1) using a GTR-gamma nucleotide substitution model. The tree with the highest log-likelihood was taken and time-calibrated using PLL-DPPDiv, which estimates divergence times using a Dirichlet Process Prior (Heath et al., 2012). DPPDiv requires a rooted phylogeny with the nodes in the correct order; however, RAxML by default produces an unrooted tree. Therefore, we included an archaeal sequence in our 16S alignment (Methanospirillum hungatei, RefSeq accession NR_074177) and used this as an outgroup in our RAxML run. This gives a tree rooted at the outgroup, which we checked for correct topology using TimeTree (Kumar et al., 2017) as a reference. We derived calibration nodes from TimeTree (Kumar et al., 2017) and performed two DPPDiv runs for 1 million generations each, sampling from the posterior distribution every 100 generations. We ensured that the two runs had converged by verifying that each parameter had an effective sample size above 200 and a potential scale reduction factor below 1.1. We summarized the output of DPPDiv into a single tree using the TreeAnnotator program implemented in BEAST (Bouckaert et al., 2019). We then dropped the outgroup tip to give a time-calibrated phylogeny of our bacterial 16S sequences only, which was used for further analysis. Details of calibration nodes used are given in Table 1. Table 1 Details of time tree calibration nodes. We constrained the time calibration of our RAxML tree based on estimated divergence times from TimeTree (Kumar et al., 2017). Taxa ATaxa BMin divergence time (MYA)Max divergence time (MYA)BacteriaArchaea4290-PseudomonasBacillus31003254PseudomonasLabrys10533135PseudomonasCollimonas10533135CollimonasVariovorax1271-ArthrobacterStreptomyces14201870ArthrobacterBacillus31003254BacillusBrevibacillus17342398 To test whether there was evidence of evolution of Topt, we calculated Pagel’s λ (Pagel, 1999), which quantifies the strength of phylogenetic signal – the degree to which shared evolutionary history has driven trait distributions at the tips of a phylogenetic tree. λ=0 implies no phylogenetic signal, that is, the signal expected if variation in trait values is independent of the phylogeny. λ=1 implies strong phylogenetic signal, that is, that the trait has evolved gradually along the phylogenetic tree (approximated as Brownian motion [BM]). Intermediate values (0<λ<1) imply deviation from the BM model, and may be observed for different reasons, such as constrained trait evolution due to stabilizing selection, and variation in evolutionary rate over time (e.g., due to episodes of rapid niche adaptation). Pagel’s λ requires that the trait be normally distributed. However, Topt in our dataset has a right-skewed distribution. Therefore, to test phylogenetic heritability we calculated λ for log⁡(Topt). Blomberg’s K is another metric that is also widely used to infer phylogenetic heritability (Blomberg et al., 2003; Münkemüller et al., 2012). Blomberg’s K calculates the phylogenetic signal strength as the ratio of the mean squared error of the tip data and the mean squared error of the variance–covariance matrix of the given phylogeny under the assumption of BM (Münkemüller et al., 2012). K=1 indicates taxa resembling each other as closely as would be expected under a BM model, K<1 indicates less phylogenetic signal than expected under BM, and K>1 indicates more phylogenetic signal than expected and thus a substantial degree of trait conservatism (Blomberg et al., 2003). Under a BM model of trait evolution, Pagel’s λ is expected to perform better than K, which may itself be better utilized for simulation studies (Münkemüller et al., 2012). Previous work suggests that Tpk is likely to evolve in a BM manner in prokaryotes (Kontopoulos et al., 2020a), making λ a more appropriate metric for these data than K. Furthermore, λ is potentially more robust to incompletely resolved phylogenies and is therefore likely to provide a better measure than K for ecological data in incomplete phylogenies (Molina-Venegas and Rodríguez, 2017). Therefore, we use λ as likely the more appropriate metric for our data; however, for the sake of completeness, we also test for phylogenetic heritability using K. We mapped the evolution of Topt onto our phylogeny using maximum likelihood to estimate the ancestral values at each internal node, assuming a BM model for trait evolution (an appropriate model, given