Enzyme adaptation to habitat thermal legacy shapes the thermal plasticity of marine microbiomes

Extraction of total active proteomes from sediment samples

We sampled 14 sediments along the coastlines of the Irish Sea, the Mediterranean Sea, and the Red Sea (from 16°N to 53°N), applying uniform sampling and storage procedures. Location details and sediment temperature fluctuations are summarized in Supplementary Table S1. We collected sediments (5 Kg) in triplicate and extracted the total proteins using a well-established microbial detachment procedure⁶⁷, with some modifications. We mixed 100 g of sediment with 300 ml of sterilized saline solution (5 mM sodium pyrophosphate and 35 g L⁻¹ of NaCl) containing 150 mg L⁻¹ of Tween 80 (from Merck Life Science S.L.U., Madrid, Spain) in an ice water bath. After re-suspension, samples were kept in a water bath ultra-sonicator (Bandelin SONOREX, Berlin, Germany) on ice and sonicated (60 W output) for 120 min. We repeated this procedure twice, with an ice water bath incubation of 60 min between each cycle. We then centrifuged the samples at 500 g for 15 min at 4 °C to remove the sediments in a centrifuge 5810 R (Eppendorf AG, Hamburg, Germany). Supernatants were carefully transferred to a new tube, minimizing disruption of the sediments, and the resulting supernatants were centrifuged at 13,000 g for 15 min at 4 °C to produce microbial cell pellets. We used the resulting cell mix to extract the total protein by mixing the cells with 1.2 ml BugBuster® Protein Extraction Reagent (Novagen, Darmstadt, Germany) for 30 min with shaking (250 rpm). Subsequently, samples were disrupted by sonication using a pin Sonicator® 3000 (Misonix, New Highway Farmingdale, NY, USA) for a total time of 2 min (10 watts) on ice (4 cycles × 0.5 min with 1 min ice-cooling between each cycle). Extracts were centrifuged for 10 min at 12,000 g at 4 °C to separate cellular debris and intact cells. Supernatants were carefully aspirated (to avoid disturbing the pellet), transferred to new tubes, and stored at –80 °C until use. The protein solution was filtered at 15 °C for 7 h using Vivaspin filters (Sartorius, Goettingen, Germany) with a molecular weight (MW) cut-off of 3,000 Da to concentrate the proteins up to a final concentration of 10 mg ml⁻¹, according to the Bradford Protein Assay (Bio-Rad Laboratories, S.A., Madrid, Spain)⁶⁸. The average total amount of proteins extracted per each 100 g of sediment was 612 µg (interquartile range, 31 µg, see details in Supplementary Fig. S2). In all cases, extensive dialysis of protein solutions against 40 mM 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES) buffer was performed using a Pur-A-LyzerTM Maxi 1200 dialysis kit (Merck Life Science S.L.U., Madrid, Spain)⁶⁹, and active proteins stored at a concentration of 10 mg ml⁻¹at –86 °C until use. As reported previously⁷⁰, 2DE was performed using GE Healthcare reagents and equipment, 11 cm IPG strips in the pH range of 3–10 and molecular weight ranging from 10 to 250 kDa (Precision Plus Protein Dual Color Standards #1610374, Bio-Rad Laboratories, S.A., Madrid, Spain). The 2-DE was performed using a validated pooling strategy⁷¹, in which proteins extracted from three independent biological replicates (i.e., sediments) were mixed in equal amounts and a total of 150 µg of protein were further loaded per gel. Staining was performed with SYPRO Ruby Protein Gel Stain (Invitrogen, Waltham, MA, USA). The two-dimensional SDS-PAGE (12% acrylamide) gels of extracted proteins are reported in Supplementary Fig. S2 (original gels in Source Data). The same protocol was applied to extract and analyse by SDS-PAGE the total active proteins extracted from sediment samples with different temperature variability levels (HTV, ITV, and LTV) collected in the Red Sea (Supplementary Table S4). The total amount of protein extracted per each 100 g of sediment is given in Supplementary Table S8. Coomassie-stained one-dimension SDS-PAGE (1-DE) gels of extracted proteins are shown in Supplementary Fig. S9 (original gel in Source Data).

Source, expression and purification of esterases and EXDOs from a wide geographical range

We recovered 83 enzymes (78 esterases and 5 EXDO) from microbial communities inhabiting marine sediments across ten distinct locations from the latitudinal transect described above: Ancona harbour (Anc), Priolo Gargallo (Pri), Gulf of Genoa, Messina harbour (Mes), Milazo harbour (Mil), Mar Chica lagoon (MCh), Bizerte lagoon (Biz), El-Max site (ElMax), Gulf of Aqaba (Aq), and Menai Strait (MS); further details are provided in Supplementary Data S3. Sources of the enzymes were the corresponding shotgun metagenomes (see Supplementary Table S3) and the metagenome clone libraries generated from the extracted DNA⁷¹. The sediment sample from the Gulf of Genoa was not used for activity tests and metaproteome analysis because no raw sample material was available; however, because of the possibility to access its shotgun metagenome (see Supplementary Table S3) and a metagenome clone library⁷², we used the sample for screening esterases to incorporate an additional latitude in our transect. In the case of Menai Strait (Irish Sea), five additional esterases were retrieved from a metagenome obtained from enriched cultures prepared with samples collected on 22^nd June 2019 from Menai Strait (School of Ocean Sciences, Bangor University, St. George’s Pier, Menai Bridge, N53°13′31.3″; W4°09′33.3”). The water temperature was 14 °C and the salinity was 32 p.s.u. Two enrichment cultures were set up at 20 °C: (i) SW: seawater enrichment with 0.1% lignin; the enrichment was set up using 50 ml of the sample as inoculum with the addition of 0.1% lignin (Sigma-Aldrich, Gillingham, United Kingdom) (w/v); (ii) AW: algal surface wash-off in seawater, enriched with 0.1% lignin; the enrichment was set up using 50 ml of surface wash-off after mixing of ca. 10 g of Fucus (brown algae) in the seawater and removal of plant tissue, 0.1% lignin (w/v) was added. After 92 days of incubation, 5 ml of each enrichment cultures were transferred into the new flask containing 45 ml autoclaved and filtered seawater with 0.1% lignin. This procedure was repeated on days 185 and 260, and the incubation was stopped on day 365. The DNA was extracted using 12 months using MetaGnome extraction kit (EpiCentre, Biotechnologies, Madison, WI, USA), sequenced on Illumina MiSeq™ platform (Illumina Inc., San Diego, CA, USA) using paired-end 250 bp reads at the Centre for Environmental Biotechnology (Bangor, UK), and sequencing reads were processed and analysed as described previously⁷³.

The screening, cloning and activity of a subset of 35 identified esterases have been reported previously⁷². The remaining 48 enzymes are reported for the first time in this study and were identified using naive and in silico metagenomic approaches, as detailed below. The environmental site from which each enzyme originated and the method employed for its identification are detailed in Supplementary Data S3. For naive screens addressing the recovery of new sequences encoding esterases and EXDO, the large-insert pCCFOS1 fosmid libraries made using the corresponding DNA samples, the CopyControl Fosmid Library Kit (Epicentre Biotechnologies, Madison, WI, USA) and the Escherichia coli EPI300-T1R strain were used. The nucleic acid extraction, construction and the functional screens of such libraries have been previously described⁷². In brief, fosmid clones were plated onto large (22.5 × 22.5 cm) Petri plates with Luria Bertani (LB) agar containing chloramphenicol (12.5 µg ml⁻¹) and induction solution (Epicentre Biotechnologies; WI, USA), at a quantity recommended by the supplier to induce a high fosmid copy number. Clones were scored by the ability to hydrolyze α-naphthyl acetate and tributyrin (for esterase activity), and catechol (for EXDO activity)^72,74. Positive clones presumed to contain esterases and EXDOs were selected, and their DNA inserts were sequenced using a MiSeq Sequencing System (Illumina, San Diego, USA) with a 2 × 150-bp sequencing v2 kit at Lifesequencing S.L. (Valencia, Spain). After sequencing, the reads were quality-filtered and assembled to generate nonredundant meta-sequences, and genes were predicted and annotated via BLASTP and the PSI-BLAST tool⁷². For in silico screens, addressing the recovery of new sequences encoding esterases, the predicted protein-coding genes, obtained after the sequencing of DNA material from resident microbial communities in each of the samples, were used. The meta-sequences are available from the National Center for Biotechnology Information (NCBI) nonredundant public database (accession numbers reported in Supplementary Data S3). Protein-coding genes identified from the DNA inserts of positive clones (naive screen) or from the meta-sequences were screened for enzymes of interest using the Blastp algorithm via the DIAMOND v2.0.9 program with default parameters (percentage of identity ≥60%; alignment length ≥70; e-value ≤1e⁻⁵)²⁹, against the Lipase Engineering sequence databases (to screen for esterases) and AromaDeg database (for EXDO)⁷⁴. Since the collection of sediments across locations experiencing different MATs was limited by our sampling capacity, to expand our range of exploration at a global scale and to validate our dataset, we added our single enzyme analysis to the seawater metagenomes retrieved from the Tara Ocean Expedition database (accession number in Supplementary Data S4). Due to the volume of sequences generated, this database provides access to a large number of enzymes, including those studied here through homology search. Esterases were selected as target sequences, and the following pipeline was used. First, we selected a sequence encoding an esterase reported as one of the most substrate-ambiguous esterases out of 145 tested (EH₁, Protein Data Bank acc. nr. 5JD4) and well-distributed in the marine environment⁷². Second, we performed a homology search of this sequence against the Tara Ocean metagenome²¹ to retrieve similar sequences, using the Blastp algorithm via the DIAMOND v2.0.9 program³⁰ (e-value <e⁻¹⁰). A total of 150 sequences encoding presumptive such enzymes from 56 different locations of the Tara Ocean Expedition were selected (Supplementary Data S4).

Once identified, the sequences encoding the wild-type enzymes—here identified and reported for the first time from all the geographically distinct locations (including the ones from the Tara Ocean Expedition)—were used as templates for gene synthesis. Genes were codon-optimized to maximize expression in E. coli. Genes were flanked by BamHI and HindIII (stop codon) restriction sites and inserted in a pET-45b(+) expression vector with an ampicillin selection marker (GenScript Biotech, EG Rijswijk, Netherlands). This plasmid, which was introduced into E. coli BL21(DE3), supports the expression of N-terminal histidine (His) fusion proteins with the final amino acid sequences of all synthetic proteins being MAHHHHHHVGTGSNDDDDKSPDP-X, where X corresponds to the original sequence of the target enzyme (Supplementary Data S3 and S4). In all cases, the soluble His-tagged proteins were produced and purified at 4 °C after binding to a Ni-NTA His-Bind resin (Merck Life Science S.L.U., Madrid, Spain), as described previously^72,74. Purity was assessed as > 98% using SDS-PAGE analysis in a Mini PROTEAN electrophoresis system (Bio-Rad Laboratories, S.A., Madrid, Spain). Purified protein was stored at –86 °C until use at a concentration of 10 mg ml⁻¹ in 40 mM HEPES buffer (pH 7.0). A total of approximately 5–40 mg of total purified recombinant protein was obtained from 1 L of culture. Supplementary Fig. S1 illustrates a schematic representation of the pipeline implemented in this work to investigate enzyme activities in a large set of marine samples, starting from samples collected (sediments) and available metagenomes.

Enzyme activity assessments

All substrates used for activity tests were of the highest purity and, if not indicated otherwise, were obtained from Merck Life Science S.L.U. (Madrid, Spain): 4-nitrophenyl-propionate (ref. MFCD00024664), 4-nitrophenyl phosphate (ref. 487663), 4-nitrophenyl β-D-galactose (ref. N1252), bis(p-nitrophenyl) phosphate (ref. 123943), benzaldehyde (ref. B1334), 2-(4-nitrophenyl)ethan-1-amine (ref. 184802-5G), pyridoxal phosphate (ref. P9255), acetophenone (ref. A10701), NADPH (ref. N5130) and catechol (ref. PHL82372). We directly tested total protein extracts for esterase, phosphatase, beta-galactosidase, and nuclease activity using 4-nitrophenyl-propionate, 4-nitrophenyl phosphate, 4-nitrophenyl β-D-galactose, and bis(p-nitrophenyl) phosphate, respectively, by following the production of 4-nitrophenol at 348 nm (extinction coefficient [ε], 4147 M⁻¹ cm⁻¹), as previously described⁶⁹. For determination: [total protein]: 5 μg ml⁻¹; [substrate]: 0.8 mM; reaction volume: 200 μl; T: 4–85 °C; and pH: 8.0 (50 mM Tris-HCl buffer). The hydrolysis of 4-nitrophenyl-propionate was used to determine, under these standard conditions, the effects of temperature on the purified esterase. Transaminase activity was determined using benzaldehyde as amine acceptor, 2-(4-nitrophenyl)ethan-1-amine as amine donor, and pyridoxal phosphate as a cofactor, by following the production of a colour amine at 600 nm (extinction coefficient, 537 M⁻¹ cm⁻¹), as previously described⁷⁵. For determination, [total protein]: 5 μg ml⁻¹; [substrates]: 25 mM; [pyridoxal phosphate]: 1 mM; reaction volume: 200 μL; T: 4-85 °C; and pH: 8.0 (50 mM Tris-HCl buffer). Aldo-keto reductase activity was determined using acetophenone as a substrate and NADPH as a cofactor, by following the consumption of NADPH at 340 nm (extinction coefficient, 6220 M⁻¹ cm⁻¹), as described⁷⁶. For determination, [total protein]: 5 μg ml⁻¹; [substrate]: 1 mM; [cofactor]: 1 mM; reaction volume: 200 μL; T: 4–85 °C; and pH: 8.0 (50 mM Tris-HCl buffer). We determined EXDO activity using catechol as substrate, by following the increase of absorbance at 375 nm of the ring fission products (extinction coefficient, 36000 M⁻¹ cm⁻¹), as previously described⁷⁴. For determination, [protein]: 5 μg ml⁻¹; [catechol]: 0.5 mM; reaction volume: 200 μL; T: 4–85 °C; and pH: 8.0 (50 mM Tris-HCl buffer). The hydrolysis of catechol was used to determine, under these standard conditions, the effects of temperature on the purified EXDOs. All measurements were performed in 96-well plates (ref. 655801, Greiner Bio-One GmbH, Kremsmünster, Austria), in biological triplicates over 180 min in a Synergy HT Multi-Mode Microplate Reader (Biotek Instruments, Winooski, VT, USA) in continuous mode (measurements every 30 s) and determining the absorbance per minute from the slopes generated and applying the formula (1). All values were corrected for nonenzymatic transformation.

Shotgun proteomics

Proteomics was performed by using total active proteins (extracted as above), which were then subjected to protein precipitation, protein digestion and Liquid Chromatography-Electrospray Ionization Tandem Mass Spectrometric (LC-ESI-MS/MS) analysis, as previously described⁷⁷. High-quality reference metagenomes corresponding to each sample (BioProject number in Supplementary Table S3) were used for protein calling, with a threshold of only one identified peptide per protein identification because False Discovery Rates (FDR) controlled experiments counter-intuitively suffer from the two-peptide rule. The confidence interval for protein identification was set to ≥95% (p < 0.05), and only peptides with an individual ion score above 20 were considered correctly identified. All protocols and experimental details, including those for mass spectrometry and measures of quality and fidelity of the datasets, have been previously described⁷⁷. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE⁷⁸ partner repository with the dataset identifier PXD039714 and 10.6019/PXD039714.

Thermal denaturation assessments of proteins through circular dichroism

Spectra were acquired between 190 and 270 nm with a Jasco J-720 spectropolarimeter equipped with a Peltier temperature controller, employing a 0.1 mm cell at 25 °C. Spectra were analysed, and T_d values were determined at 220 nm between 10 °C and 85 °C at a rate of 30 °C per hour in 50 mM Britton and Robinson buffer at pH 8.5. A protein concentration of 1.0 mg ml⁻¹ was used. T_d (and standard deviation of the linear fit) was calculated by fitting the ellipticity (millidegrees; θ) at 220 nm at each of the different temperatures using a 5-parameter sigmoid fit with SigmaPlot 13.0^79,80. We used a generalized additive model (GAM)⁸¹ to analyse T_d (our response variable) at each temperature (continuous explanatory variable) in the hot and cold seasons (in August and December, respectively) and the three levels of temperature variability (HTV, LTV, ITV) as our two categorical explanatory variables.

Enzyme structure prediction, ensemble generation and rigidity analysis

We applied the AlphaFold2-based workflow of ColabFold^82,83, which is available at https://github.com/sokrypton/ColabFold (accessed 22.02.2022), to generate 3D structural models of esterases. A single model was generated for each esterase with ten prediction cycles (–num_recycles) and structurally refined by running a relaxation with AMBER (–amber). For subsequent analyses, only esterases with a sufficient 3D structural model quality, with a sequence length <1000 residues, and without cofactors were considered. To test whether the catalytically active residues (CARs) of the 3D structures are accessible for substrates, we used the CAVER 3.0.3 PyMOL Plugin^84,85. Therefore, CARs were identified based on the minimal summed distances between all triplets of serine oxygen atoms, histidine nitrogen atoms in epsilon position, and carbon atoms of the carboxylic acid group from either glutamate or aspartate residues. Starting points for the computations were defined based on the Cartesian coordinates of the centre of mass (COM) of each CAR. Default values were used for the probe radius (0.9 Å), shell radius (3.0 Å) and shell depth (4.0 Å). We verified that CARs in all models are accessible for substrates, i.e., that all models are in an open conformation: CARs are either located on the protein surface or are buried and connected with the surface by tunnels.

The esterase structures were pre-processed with pdb4amber, which is part of AmberTools21, and hydrogen atoms were added using the Reduce program⁸⁶. The prepared esterases were solvated in a truncated octahedron of TIP3P water⁸⁷, leaving at least 20 Å between the esterase structure and the edges of the solvent box, using the LeaP program of AmberTools21. All systems were neutralized by adding Na⁺ or Cl⁻ ions as needed. We used the Amber ff14SB force field⁸⁸ to parametrize the protein. Ion parameters were taken from Joung and Cheatham⁸⁹. Structural ensembles of esterases were generated by all-atom molecular dynamics (MD) simulations, with five replicas at 500 ns, yielding 2.5 μs of cumulative simulation time per esterase. Minimization steps, thermalization, and production simulations were carried out using the GPU-accelerated CUDA version of PMEMD^90,91 from the Amber21 suite of programs⁹². The systems were heated to 298 K and the pressure was adapted in NPT simulations to obtain a density of 1 g cm⁻³. During thermalization and density adaptation, we kept the solute fixed by positional restraints of 1 kcal mol⁻¹ Å⁻², which were gradually removed over five steps in short subsequent NVT simulations. Afterwards, five NVT production simulations of 500 ns length were performed using unbiased MD simulations. During these simulations, we set the time step to integrate Newton’s equation of motion to 4 fs, following the hydrogen mass repartitioning strategy⁹³. Coordinates were saved every 200 ps yielding 2500 conformations for each production run that were considered for subsequent analyses. The Amber21 software for MD simulations is available at http://ambermd.org/.

Rigidity analyses were performed following Nutschel et al.⁴¹ using the Constraint Network Analysis (CNA) software package (version 3.0)^33,34,94,95, available at https://cpclab.uni-duesseldorf.de/index.php/Software. In detail, we applied CNA on ensembles of network topologies generated from conformational ensembles obtained from MD simulations generated at room temperature, for which force field parameters had been optimized. Average stability characteristics were calculated by constraint counting on each topology in the ensemble. CNA functions as a front- and back-end to the graph theory-based software Floppy Inclusions and Rigid Substructure Topography (FIRST)⁹⁶. The application of CNA to biomolecules aims to identify their rigid cluster and flexible region composition, which can aid in understanding the biomolecular structure, stability and function. As the mechanical heterogeneity of biomolecular structures is intimately linked to their diverse biological functions, biomolecules generally show a hierarchy of rigid and flexible regions³⁷. To monitor this hierarchy, CNA performs thermal unfolding simulations by consecutively removing noncovalent constraints (hydrogen bonds and salt bridges) from a network in the order of their increasing strength. Therefore, a hydrogen bond energy E_HB is computed from an empirical energy function⁹⁷. For a given network state σ = f(T), hydrogen bonds (including salt bridges) with an energy E_HB > E_cut(σ) are removed from the network at temperature T. In the present study, thermal unfolding simulations were carried out by reducing E_cut from −0.1 kcal mol⁻¹ to –6.0 kcal mol⁻¹ with a step size of 0.1 kcal mol⁻¹. E_cut can be converted to a temperature T using the linear equation introduced by Radestock et al.³³ as reported in equation (2), where the range of E_cut is equivalent to increasing temperature from 302 K to 420 K with a step size of 2 K.

The number of hydrophobic tethers was kept constant during the thermal unfolding simulations⁹⁸. From these simulations, CNA computes a set of indices to quantify biologically relevant characteristics of protein stability at a global and local scale^34,99. Here, we used the cluster configuration entropy H_type2, a measure of the global structural stability, to predict the phase transition temperature T_p (for details on H_type2^34,37,99). At T_p, the protein switches from a rigid (structurally stable) to a floppy (unfolded) state and the largest rigid cluster no longer dominates the whole protein network. If the largest rigid cluster dominates the whole protein network, H_type2 is low because of the limited number of possible ways to configure a system with a very large cluster. When the largest rigid cluster starts to decay or stops to dominate the network, H_type2 jumps. At this stage, the network is in a partially flexible state with many ways to configure a system consisting of many small clusters. The percolation behaviour of protein networks is usually complex, and multiple phase transitions can be observed^{33,35,37,38,39,40,41,100}. To identify T_p, a double sigmoid fit was applied to an H_type2 versus T(E_cut) curve as performed previously^{33,35,37,38,39,40,41,100}. In general, T_p was taken as the T value associated with the largest slope of the fit, except for esterases with T_d > 50 °C for which the second phase transition was chosen to focus on the decomposition of the core. It is important to note that applying CNA to MD simulations at room temperature may lead to an evening out of T_p values for esterases that transition around this temperature, i.e., systems with a T_p at or below room temperature might all be influenced similarly by loosening their bonding network. By contrast, systems with a transition temperature at or above room temperature would still be discriminated against. The data generated in this study for analyzing T_p values have been deposited at researchdata.hhu.de under accession code DOI: 10.25838/d5p-42¹⁰¹ [https://doi.org/10.25838/d5p-42].

Relationship of temperature-induced changes in enzyme

Relationship between MAT and enzyme response to temperature (i.e., T_opt, T_d and T_p) were evaluated by performing linear regression in R. In the case of enzymes retrieved from the Tara ocean dataset we calculated first the break point (flexus) using the package segmented in R¹⁰² and then we computed separately the linear model describing the two linear regressions before and after the breakpoint. To evaluate the possible relation between enzyme thermal response and other environmental parameters, salinity and pH data were retrieved from Bio-ORACLE⁵² using GPS coordinates of each location.

Environmental characterization and sediment collection from different temperature variability levels in the Red Sea

We recorded the temperatures of surface sediments from March 2015 to September 2016 along the coast of the Red Sea using HOBO data loggers (Onset, USA) in nine stations located at 3, 25, and 50 m depth. Details on the location, depth and temperature fluctuations of the studied sediments are reported in Supplementary Table S4 and Source Data. We first assess the differences in the homogeneity of the temperature variance in the three types of sediments to evaluate the magnitude of thermal variation and then we test the difference among their MATs using a non-parametric ANOVA (Dunnett’s multiple comparisons tests). We identified three different levels of temperature variability (Fig. 3a–c; Supplementary Table S5): high, intermediate, and low thermal variability (HTV, ITV, and LTV, respectively), where sediments experienced temperature variations of 12.8 °C, 8.8 °C, and 6.7 °C, respectively. From each station, we sampled 200 g of surface sediment (0–5 cm depth) in triplicate in August and December 2015 with a Van der Venn grab (1 dm³) equipped with a MicroCat 250 Seabird CTD (Conductivity, Temperature, Depth), which was assembled on board the research vessel R/V Explorer (KAUST). During sampling, we measured the temperature of the sediments and the water layer covering the sediments using a digital thermometer and the CTD, respectively. We conducted all sampling in compliance with the guidelines specified by KAUST and Saudi Arabian authorities.

Sediment processing for analysis of bacterial communities

From each sample (in triplicate), we immediately removed subsamples of sediment (n = 54, ~10 g) and stored them at –20 °C for molecular analysis. Separately, sediment 25 ± 1 g was transferred to 50 ml tubes and added 30 ml of filtered (0.2 µm) water from the Red Sea. The tubes were shaken at 500 rpm for one hour and then centrifuged them at 300 g for 15 min to detach the microbial cells in the sediments without affecting their vitality^103,104. The supernatant containing the extracted cells was collected in sterile tubes and was immediately used to measure microbial growth rates.

Evaluation of bacterial growth in sediments at different temperatures

We evaluated the microbial growth rate of the heterotrophic community extracted from the sediments under HTV, ITV, and LTV at 10 °C, 20 °C, 30 °C, 40 °C and 50 °C, using Marine Broth as the cultivation medium (Zobell Marine Broth 2216) supplemented with 0.1 g/L cycloheximide; a rich-medium was selected to avoid the nutrient limitation effect that can affect bacterial physiology^63,105. We inoculated 96-well plates with 200 µl of cultivation medium and 25 µl of the cell suspension extracted from the sediments. We inoculated the three biological replicates from each station and each level of temperature variability in eight wells, giving a total of 72 wells for each plate, with 24 wells used as a negative control inoculated with water. We assembled a total of three plates for each incubation temperature from August and December. Plates were spectrophotometrically measured at 3 h intervals using an optical density of 600 nm (Spectramax® M5) for 72 h. Wells with optical density <0.15 (average value of the medium inoculated with autoclaved sediment extracts) were assigned as ‘no-growth’ at a given temperature. Growth (expressed as OD₆₀₀) was normalised in terms of the OD₆₀₀ of the initial inoculum. The growth rate was calculated as the change in the number of cells in a culture per unit of time (h)¹⁰⁶. A GAM⁸¹ was applied to evaluate the effect of temperature, temperature variability level and seasonality on the continuous response variables of bacterial growth and growth rate.

Total DNA extraction, Illumina sequencing, and taxonomic analysis of bacterial 16S rRNA gene sequences

The total DNA was extracted from 0.4 ± 0.1 g of sediment using a DNeasy PowerSoil Pro Kit® (Qiagen) and from the final enriched heterotrophic bacterial fraction obtained by the cultivation approach using a DNeasy UltraClean Microbial Kit (Qiagen). PCR amplification of the V3–V4 hypervariable regions of the 16S rRNA gene on DNA in the sediment samples was performed using the universal primers 341f and 785r¹⁰⁷. We constructed libraries with the 96 Nextera XT Index Kit (Illumina) following the manufacturer’s instructions and sequenced DNA using the Illumina MiSeq® platform with paired-end sequencing at the Bioscience Core Lab at KAUST. Raw reads were deposited in the NCBI database under the SRA accession number PRJNA508596. We assembled forward and reverse reads for each sample into paired-end reads (minimum overlap of 50 nucleotides and a maximum of one mismatch within the region) using the fastq-join algorithm, and the samples were analysed using the UPARSE v8 and QIIME v1.9 softwares¹⁰⁸. The final reads (average length of 405 bases) were clustered into operational taxonomic units (OTUs), taking 97% sequence identity as the cut-off. All samples showed a sufficient sequencing depth (Good’s coverage values >90%) for further analysis (Supplementary Tables S9 and S10). We calculated the compositional similarity matrix (Bray-Curtis of the log-transformed OTU table) with Primer 6¹⁰⁹. Using the same software, canonical analysis of principal coordinates (CAP)¹¹⁰ was used to compare the temperature variability samples (temperature variability levels: HTV, ITV, and LTV; season levels: August and December) based on the compositional similarity matrix. We applied permutational multivariate analyses of variance to the matrix (PERMANOVA; main and multiple comparison tests). We tested the occurrence of thermal-decay patterns in sediments with different temperature variability levels using linear regression (Prism 9.2 software, La Jolla California USA, www.graphpad.com) between the bacterial community similarities (Bray-Curtis) and the temperature differences among sediments (∆T°C) at the time of sampling. We calculated alphadiversity indices (richness and evenness) using the paleontological statistics (PAST) software, and their correlation with temperature was modelled using linear regression in Prism 9.2. Spearman correlation among temperature and relative abundance of OTUs within each sediment sample was evaluated; OTUs were classified based on their positive (enriched) and negative (depleted) correlation with sediment temperature.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com