Thermal adaptation of soil microbial growth traits in response to chronic warming

All organisms were isolated from soils collected from the Harvard Forest long-term warming study, located at the Harvard Forest LTER site in Petersham, MA (21). The site is a mixed hardwood forest with paper and black birch (Betula papyrifera and lenta), red maple (Acer rubrum), black and red oak (Quercus velutina and rubra), and American beech (Fagus grandifolia) dominant tree species. Soils are coarse-loamy inceptisols. Eighteen 6 × 6 m² plots were randomly assigned one of three treatments: (i) plots with buried electrical cables, heating soils 5°C above ambient temperature throughout the year; (ii) disturbance control plots with the same as set up as the heated plots, but without electrical power; and (iii) undisturbed control plots. Soils are heated 5°C above ambient temperature by way of electrical cables buried 10 cm below the soil surface. This temperature was chosen as a worst-case scenario rise in soil temperatures by the year 2100 (1).

We selected 23 strains of Alphaproteobacteria (Table 1) from our lab culture collection originating from either the heated (n = 8) or control plots (n = 15). Bacteria were isolated from soils using several cultivation methods (Table S1) and cryopreserved at −80°C. Isolates were grown on 10% Tryptic Soy Agar until we could identify distinct colony morphology. We genotyped isolates by sequencing their full-length 16S ribosomal RNA. We extracted genomic DNA using CTAB-lysozyme extraction protocol (22). 16S rRNA was amplified on an Eppendorf AG 22331 Hamburg using the 27F and 1492R primers. We used a 25-µL final reaction volume with 0.125 µL Invitrogen Taq, 10 µL MgCl₂ 10× PCR buffer, 0.75 µL 50 mM MgCl₂, 1 µL of each primer, 2 µL of dNTP mix, and 1 µL of template for amplification reactions. We performed PCR amplifications using 35 cycles of 94°C (45 s), 50°C (30 s), and 72°C (120 s), followed by a final extension of 72°C (10 min). We used agarose gel electrophoresis to verify amplifications. DNA purification and Sanger sequencing were performed by Genewiz at Azenta Life Sciences (15).

To extract genomic DNA for genome sequencing, strains were grown in 10% Tryptic Soy Broth, and pellets were extracted using the DNeasy Blood and Tissue kit (Qiagen). We let cultures grow until the late exponential phase (OD₆₀₀ nm 0.6–0.8). One day before extractions, we added 100 µL of 10% glycine to culture tubes for a final concentration of 1%; this helped to prevent cells from adhering to one another and forming clumps as cultures reached the late exponential phase. DNA was eluted in TE buffer, quantified by Qubit, and transferred to the freezer for long-term storage.

Genomes were sequenced by the United States Department of Energy’s Joint Genome Institute (JGI), the University of Massachusetts Medical Center, or in-house using an Oxford Nanopore Technologies (ONT) MinION (Table S1). Illumina sequencing technology was performed at JGI and UMass Medical Center according to standard operating procedures (23). Long-read ONT libraries were prepared with the Ligation Sequencing Kit SQK-LSK-109 and samples were multiplexed using the Native Barcoding Expansion Kit EXP-NBD104 (Oxford Nanopore Technologies, UK). The Oxford Nanopore Native Barcoding Protocol (Oxford Nanopore Technologies) was followed, and 6–8 strains were multiplexed together in a run. Genomes were resequenced until at least 100X coverage was reached. The Covaris g-TUBE shearing step was skipped to target long fragment DNA. Starting with 1 µg of DNA per strain, samples were repaired and end-prepped using the NEBNext FFPE DNA Repair Mix and NEBNext Ultra II End Repair/dA-Tailing kits (New England Biolabs, USA). DNA was cleaned using Ampure XP Beads (Beckman Coulter, USA). Samples were ligated to individual barcodes, and then 150 ng of each sample pooled together, for a final library of 700–1,000 ng. Adapters were ligated to the sample with Blunt/TA ligase (New England Biolabs). The long fragment buffer provided in the sequencing kit was used in an extended 10 min incubation at 37℃ to enrich for high-molecular weight DNA. The flow cell was primed using the Flow Cell Priming Kit (Oxford Nanopore Technologies), and 15 fmol of the library was mixed with the sequencing buffer and loading beads and then loaded through the Spot-On port of the flow cell (24).

For genomes sequenced using ONT, sequence runs were initially basecalled using the high-accuracy base calling (HAC) algorithm with Guppy (25). These fast5 files were concatenated into one file, then reads were subsampled based on read quality using Filtlong (26). The filtered reads were de novo assembled using Flye (27). A consensus assembly was generated using Racon (28), and final polishing was performed with Medaka (29). These assemblies were checked for quality using Quast (30) and CheckM (31).

We measured bacterial growth over time using absorbance measurements for liquid cultures spanning temperatures from 22°C to 37°C along 2–3°C increments. Because organisms tend to live at or below their optimal growth temperature (32), we chose 22–37°C to capture the temperature growth range of mesophiles, the temperature range during the growing season at the Harvard Forest (33), and to accommodate instrumental limitations. Absorbance was measured in 96-well plates by optical density at 600 nm (OD₆₀₀ nm) as a measure of cell abundance. Each well was filled with 240 µL of 10% Tryptic Soy Broth. Colonies grown on petri plates were resuspended in 500 µL of phosphate-buffered saline and 10 µL resuspension was inoculated into each well. For temperatures at or above 30°C, we pipetted 2–3 mL of a 0.05% solution of Triton X-100 in 20% ethanol on the plate lid to prevent condensation. Each plate accommodated eight isolates with 11 replicates each and eight negative controls according to a randomized plate format. Bacterial growth was measured using a SpectraMax M2 plate reader (Molecular Devices, CA, USA) at OD₆₀₀ nm. Growth curves lasted 72–99 h or until microbes entered death phase.

We fitted the Gompertz growth curve on data for OD₆₀₀ nm over time to calculate growth rate using the R package Growthcurver (34) (Fig. 1A). A Gompertz growth curve is an established time course model that parametrizes bacterial growth over time as a sigmoidal function (17). We estimated the intrinsic growth rate from the fitted Gompertz model. This was repeated for each replicate at each temperature.

We fitted temperature response curves to estimate the microbial growth traits. Temperature sensitivity of growth is an estimate of how growth rate changes with increasing temperature, but it can be estimated using different parameters depending on the model applied (Table 1). Optimum growth temperature is the temperature at which growth rate is the greatest. Maximum growth temperature is the estimated highest temperature at which microbial growth occurs. We modeled the relationship between growth rate and temperature for each isolate using the Ratkowsky 1983 model (Fig. 1B):

T _min is the minimum permissible temperature for growth (°C), T _max is the maximum permissible growth temperature (°C), and c is an empirical parameter required to model data above the optimum temperature (°C⁻¹). Temperature sensitivity of growth was quantified by the Ratkowsky parameter b (OD₆₀₀nm*min^−0.5/°C), which is the regression coefficient for square root of growth rate on temperature (18, 35). Maximum growth T _max (°C) was extracted as a parameter from the fitted model. Optimum growth Topt (°C) was estimated from the fitted model (36). Model fitting was performed using the R package rTPC (36).

The temperature optima (Topt) and inflection point (Tinf) of each bacterial growth curve were estimated using a modified version of MMRT (Fig. 1C):

where k is the bacterial growth rate (OD₆₀₀nm), k _B is Boltzmann’s constant, T is the temperature (K), h is Planck’s constant, R is the universal gas constant, $\Delta H_{T0}^{‡}$ (superscript denotes transition state) is the change in enthalpy (J mol⁻¹), $\Delta S_{T0}^{‡}$ is the change in entropy (J mol⁻¹ K⁻¹), $\Delta C_P^{‡}$ is the change in heat capacity (J mol⁻¹ K⁻¹), and T ₀ is the reference temperature (set to 296.1K) (37). Temperature inflection point is the temperature at which the greatest change in growth rate occurs. The MMRT equation was modified to allow $\Delta C_P^{‡}$ to vary linearly with temperature:

where A is the slope and B is the value of $\Delta C_P^{‡}$ at T ₀ (20). We used a non-linear least-squares regression in R version 4.2.1 to fit MMRT and calculated the Topt and Tinf numerically using the first and second derivatives, respectively. We chose to use the modified version of MMRT to better capture the Topt due to asymmetries observed in the temperature response data and because $\Delta C_P^{‡}$ varies over a wide temperature range (38 – 40). We calculated residual standard errors to determine whether the modified MMRT or original MMRT more adequately fit our data (Table S2).

To evaluate the fit of the models, we calculated the residual standard error (RSE). The Ratkowsky parameter (b), Topt estimated by both Ratkowsky 1983 and MMRT, T _max estimated by Ratkowsky 1983, and Tinf estimated by MMRT were used as traits for the phylogenetic group comparison. Outliers, potentially due to irregularities in the replicate, were excluded from analysis if they were not within the same order of magnitude as the remaining points in the data set. No more than six to seven outliers were removed. Outliers were typically due to experimental noise at higher temperatures, which was determined by assessing replicate-specific growth curves.

We conducted a phylogenetic group comparison of traits (11 – 13) to test our hypotheses that Alphaproteobacteria from warmed plots have (i) less temperature-sensitive growth rates; (ii) higher optimum growth temperatures; and (iii) higher maximum growth temperatures compared to isolates from control plots. We used the R package nlme to conduct PGLS test (41, 42). A phylogenetic group comparison accounts for the lack of independence in phylogenetic hierarchical species data. A Lilliefors test for normality was used to determine whether residuals of PGLS tests were normally distributed, and Q-Q plots were made. Trait data were transformed if the Lilliefors test failed. We removed outlier data points, which were three to four orders of magnitude greater than the remaining points in the data set.

A genome-based phylogeny was constructed using the United States Department of Energy’s Systems Biology Knowledgebase (KBase) (43) (Fig. 2). Genomes were annotated by Prokka (44). The phylogeny was constructed using Insert Genome into SpeciesTree v2.2.0 (45), which creates a multiple sequence alignment based on universal genes defined by COG (Clusters of Orthologous Groups) gene families. We set the nearest public genome count to one and removed the public node in R using the package ape (46). Temperature sensitivity of growth, optimum growth temperature, and maximum growth temperature were mapped as traits on the phylogeny.

We also calculated phylogenetic signal using Pagel’s λ to quantify the tendency of the isolates to closely resemble each other based on their phylogenetic distribution of traits. Pagel’s λ is a measure of correlation between species under Brownian Motion (47) and was estimated using the R package phytools (48). P values and group means were calculated for each microbial growth trait from warmed and control soils. Soil microbial isolates were the experimental unit, and all P values less than 0.05 were considered statistically significant relationships. All statistical analyses were performed in R using RStudio (49).

There was a significant difference between optimum growth temperature (Topt) of isolates from heated versus control plots quantified by the Ratkowsky 1983 model (t(23) = 2.84, P = 0.01, n _warm = 8, n _control = 15) (Table S3). Optimum growth temperature for isolates from warmed plots (M = 30.59, SD = 1.65) was greater than those of control plots (M = 29.75, SD = 1.88) (Table 2). Residuals were normally distributed according to the Lilliefors test (P = 0.41, n _warm = 8, n _control = 15). Pagel’s λ showed that Topt was not distributed according to Brownian motion (λ = 6.61E−05, P = 1.00) (Fig. 3A).

Our results showed no significant difference in temperature sensitivity of growth quantified by the Ratkowsky parameter (t(20) = 1.79, P = 0.09, n _warm = 7, n _control = 13) (Table S3). Temperature sensitivity of isolates from warmed and control plots were 0.0049 (SD = 0.001) and 0.0047 (SD = 0.002) respectively (Table 2). Residuals were normally distributed according to the Lilliefors test (P > 0.05). Pagel’s λ showed that temperature sensitivity of growth was not distributed according to Brownian motion, suggesting a random distribution (λ = 6.61E−05, P = 1.00) (Fig. 3B).

Results of the phylogenetic least squares test for maximum growth temperature quantified by the Ratkowsky 1983 model showed no significant difference between isolates from heated and control plots (t(23) = −0.35, P > 0.05, n _warm = 8, n _control = 15) (Table S3). Maximum growth temperature for isolates from the warmed and control plots were 37.17°C (SD = 1.00) and 37.06°C (SD = 0.88), respectively (Table 2). Residuals were normally distributed according to the Lilliefors test (P = 0.45). Pagel’s λ showed that maximum growth temperature was not distributed according to Brownian motion (λ = 6.61E−05, P = 1.00) (Fig. 3C).

Results of the phylogenetic generalized least squares test for optimum growth temperature quantified by MMRT showed no significant difference between isolates from the heated and control plots (t(23) = 1.60, P = 0.12, n _warm = 8, n _control = 15) (Table S4). Optimum growth temperature for isolates from warmed plots (M = 30.26, SD = 1.82) was not significantly different than those of control plots (M = 30.74, SD = 1.72) (Table 2). Residuals were normally distributed according to the Lilliefors test (P > 0.05). Pagel’s λ showed that Topt was not distributed according to Brownian motion (λ = 0.45, P > 0.05) (Fig. 4). Residual standard errors indicated that the Ratkowsky 1983 model more adequately fit our data in comparison to the MMRT model (Table S2).

There was no significant difference in temperature inflection point (Tinf) between isolates from the heated and control soils (t(23) = 1.04, P > 0.05, n _warm = 8, n _control = 15) (Table S4). Residuals of PGLS on untransformed Tinf data failed the Lilliefors test for normality (P < 0.05). Average temperature inflection point for isolates from warmed plots was 28.00°C (SD = 2.84). Average temperature inflection point for isolates from control plots was 26.73°C (SD = 2.45) (Table 2). Residuals failed the Lilliefors test for normality following log, square root, and lambda (Box-Cox test) transformations (P < 0.05). Pagel’s λ showed that Tinf was not distributed according to Brownian motion (λ = 6.61E−05, P > 0.05) (Fig. 4). Since the residuals failed the Lilliefors test for normality, we transformed Tinf (i.e., log transformation, square root, and box-cox). However, all transformations also failed the test for normality, and suggest that results of PGLS for Tinf estimated by MMRT should be interpreted with caution.