Legacies of temperature fluctuations promote stability in marine biofilm communities

Study site

The study was performed along the coast of Calafuria (Livorno, 43° 30’N, 10°19’ E) between May and September 2018. The coast is composed of gently sloping sandstone platforms with high-shore levels (0.3–0.5 m above mean low water level) characterized by assemblages of barnacles interspersed among areas of seemingly bare rock, where the biofilm develops.

Experimental design

To assess how different levels of warming fluctuations affected the response of biofilm to subsequent extreme events, we imposed two temporal regimes of warming: 1) fixed, where we imposed pulses of temperature of the same magnitude (12 °C above ambient temperature), 2) fluctuating, in which temperature varied at each pulse application, producing a scenario with the same mean temperature as the fixed-warming scenario (12 °C above ambient temperature), but with a different variance (SD = 5). To generalize the effect of the fluctuating warming treatment (i.e., not making it contingent on one specific sequence of events), we established three different temporal warming sequences (fluct-s1, fluct-s2, fluct-s3) that had the same mean (ΔT = 12 °C) and variance (SD = 5) but different temporal profiles. Warming manipulations were performed six times at a biweekly frequency from May to August 2018. The warming session consisted of maintaining for 70 min the difference between the warming chamber and the ambient air temperature as close as possible to a nominal treatment level (e.g., +12 °C). The warming device consisted of an aluminum box equipped with an adjustable butane stove and an opening to regulate the heating (Supplementary Fig. 1b). According to previous studies^49,50, an increase of 12 °C relative to the ambient temperature was expected to reduce biofilm biomass without completely eradicating biofilm. According to study-site data, warming pulses represent moderately recurrent thermal events typical of midday hours, especially during low tide and periods of high barometric pressure (Fig. 1c). Our aim was not to strictly mimic temperature predictions from regional models but rather to simulate an increase in temporal variability, which is likely to become increasingly prevalent in the Mediterranean area due to ongoing climate change^51,52,53.

To assess the stability of biofilm, after the first series of warming treatments (phase 1 of the experiment), an extreme thermal perturbation was applied in which the biofilm was exposed to a temperature of 60 °C for 120 min twice in four days (phase 2). This was a strong, but realistic perturbation, as temperature peaks of 60 °C for 2 h in sunny conditions have been occasionally recorded at the study site^49,51 (Fig. 1c). To determine the effect of extreme events, we established four additional plots that received only the extreme temperatures and were therefore not previously disturbed (Extreme-Only). The warming device consisted of an aluminum box equipped with an adjustable butane stove and an opening to regulate the heating (Supplementary Fig. 1b). Aerial temperature was constantly measured with iButton loggers during the warming sessions inside and outside the heating chambers. Experimental plots were located 2–8 m apart and consisted of portions of substratum of 40 × 40 cm marked at their corners with rawlplugs inserted into the rock for subsequent relocation. Four experimental plots were randomly assigned to experimental conditions and interspersed across the study area. Four unmanipulated plots were used as controls. Three additional plots were established as control for artifacts (CA) to assess the potential effect of shading on biofilm during the heating sessions. CA plots were shaded with cardboard chambers of the same size of warming ones, but they were not warmed.

Data collection

Biofilm biomass was determined by means of an image-based remote sensing technique that uses chlorophyll a concentration as a proxy for biofilm biomass (μg chl a cm⁻²)²⁰. Chlorophyll a was estimated from a ratio of reflectance at near-infrared (NIR) and red bands (Ratio Vegetation Index – RVI) by means of a IR-sensitive camera (ADC)¹⁹. The exposure time of the camera was set according to the light conditions that existed at the moment of image acquisition. Because light conditions varied during sampling, a reflectance standard (30% reflective Spectralon®), with a Lambertian surface, was placed within the field of view of the camera. NIR/red ratios were related to chlorophyll a by a linear relationship, calculated based on of laboratory chlorophyll a extraction from Calafuria sandstone cores²⁰. For each photo of a plot, six subplots of 128 × 128 pixels were haphazardly taken and then processed with a Java routine on ImageJ software to obtain biofilm biomass estimates for each pixel. The spatial resolution (i.e., the area of each pixel) was 0.2 mm². Biofilm biomass was evaluated in nine times during the duration of the study, except for Extreme-Only plots, which were sampled starting from mid-June, 40 days before the imposition of temperature extremes. Notably, we ensured that subplots were spaced at a minimum distance of 30 cm from one another to minimize the likelihood of reselecting the same subplot on different sampling dates. Subplots were pooled across replicates and used exclusively for visualization purposes in Fig. 6b, while all statistical analyses were conducted using replicate-level data (plot) (see below).

Whole-metagenome sequencing and assembly of metagenome assemble genomes

We collected surface rock samples (4 × 4 cm) covered in biofilm with a hammer and chisel from three plots per treatment at each of four time points: at the beginning of the experiment, after the imposition of thermal history, after the extreme events, and at the end of the experiment. Samples were taken haphazardly, leaving a buffer area of 5 cm from the edges of neighboring samples. A temporal lag of 5 to 7 days between treatment application and sample collection allowed the biofilm to adequately respond to thermal perturbations. Rock samples were kept on ice and in dark conditions for less than two hours before being stored at −80 °C. Prior to DNA extraction, we manually removed the excess rock to increase the biofilm/rock ratio and maximize the yield of the DNA extraction. The rock samples were then pulverized with a mortar and pestle, and a 0.5 g sample was taken for DNA extraction through the Qiagen DNeasy PowerSoil Kit according to the manufacturer’s instructions. DNA concentration and quality were determined using both NanoDrop ND-1000 spectrophotometer and Qubit Fluorometric Quantitation (DS DNA High-Sensitivity Kit, ThermoFisher). DNA sequencing libraries were prepared at IGA Technology Services Srl (Udine, Italy) and then sequenced using an Illumina HiSeq2500 sequencer (Illumina, San Diego, CA) in a 2 × 150 bp paired-end mode, generating a total of 959,137,528 reads, with an average of 12,456,332 reads per sample. Reads were de-multiplexed based on Illumina indexing system. One sample (fixed warming at the initial stage of the experiment) could not be sequenced due to an extremely low DNA concentration, resulting in a total of 77 samples available for analysis.

To reconstruct the Metagenome-Assembled Genome (MAG), we used the MetaWRAP pipeline v1.3.2⁵⁴. Using the read_qc module with default parameters, raw reads were quality-trimmed, and human contamination was removed. Only high-quality reads from each sample were individually assembled into contigs through MEGAHIT (v1.1.2)⁵⁵ with a minimum contig length of 1000 bp; the quality of these assemblies was evaluated using QUAST (v5.0.2)⁵⁶. High-quality contigs were grouped into Metagenome-Assembled Genome (MAGs) using the maxbin2, metabat2, and concoct algorithms included in MetaWRAP⁵⁴. The output from these three algorithms was then integrated and refined into an optimized set of MAGs from a single assembly. As a result, 137 MAGs were generated, and their quality was assessed using CheckM (v1.0.18)⁵⁷. Finally, the relative abundance of each MAG in each sample was calculated as the length-weighted average of the contigs in each MAG and expressed as Transcripts per Million (TPM).

Taxonomic annotation of MAGs was performed using GTDB-Tk v1.3.0, which is based on GTDB ver 95⁵⁸. PhyloPhlAn 3.0 was used to create a maximum-likelihood phylogenetic tree for 136 MAGs⁵⁹. One MAG (mag_18) was excluded from the analysis due to excessive fragmentation, which led to poor alignment and the detection of fewer than 50 marker genes, making its phylogenetic placement unreliable.

To address the potential underestimation of rare taxa in MAG reconstruction and better characterize the taxonomic diversity of biofilm, we first extracted 16S rRNA gene sequences from the quality-trimmed reads using SortMeRNA v4.2⁶⁰, with the SILVA ver. 128.2⁶¹ reference database, which includes both bacterial and archaeal sequences. These reads were then assembled using SPAdes v3.15.2⁶² in metagenomics mode, tailored for the diversity and complexity of metagenomic datasets, yielding a set of contigs representing reconstructed 16S rRNA gene sequences. We aligned quality-trimmed paired-end reads from each sample back to these 16S contigs using Bowtie2 v2.4.4⁶³ with default settings. The abundance of each contig was quantified by counting the aligned reads, with read counts normalized to TPM to account for variations in contig length and sequencing depth across samples.

Functional profiling

We used a novel analytical framework (microTrait)^64,65 to map and translate microbial genomes into a hierarchical trait space that included energy acquisition, resource acquisition, stress tolerance, and life history traits that underlie different microbial strategies⁶⁶. This analysis involved two main steps: 1) genome sequences were converted into open reading frames using Prodigal⁶⁷ the resulting protein sequences were then scanned through a library of gene-level Hidden Markov Models (microTrait-HMMs, with HMMER/hmmsearch from http://hmmer.org) to produce a count table for the detected gene modules. These procedures generated trait matrices (MAG × ecological traits) at different resolutions corresponding to the levels of the functional hierarchy. MicroTrait rules mapped protein family content into traits using Boolean logic. Each protein family is represented as a Boolean variable (equals to 1 if detected, 0 otherwise) whose value was determined by the output of the corresponding microTrait-HMM. Cross-references to KEGG orthologs (KO, from https://www.genome.jp/kegg-bin/get_htext) (when available) were used to determine the cut-offs.

The results of microTrait analysis should be interpreted cautiously as they represent the potential functions of a microbial community. microTrait established strict gene-level microTraitHMM thresholds based on KEGG Orthology (KO) group genes (with an accuracy greater than 99.99%), which, combined with high levels of genome completeness, resulted in highly accurate trait assignment. Additionally, it focuses on mechanistically well-studied traits with documented genetic bases, providing more accurate functional trait characterization compared to methods like Picrust2, which relies on predicted metagenomes derived from Amplicon Sequencing Variants (ASVs).

The minimal doubling time (MDT; hours) of each MAG was calculated from codon usage bias patterns in each MAG with >10 ribosomal proteins using the R package “gRodon”^26,68.

Data analysis

Taxonomic and functional diversity

Taxonomic diversity was estimated as the equivalent number of species by calculating the first three Hill numbers using the hillR R package: MAG richness (q = 0), the exponential of Shannon’s entropy (q = 1, referring to Shannon diversity), and the inverse of Simpson’s concentration (q = 2, referring to Simpson diversity)⁶⁹. Spatial synchrony was calculated for each plot as the average population synchrony across MAGs, using Gross’ η as the synchrony metric, weighted by their relative abundance^70,71. Population variability was obtained as the weighted mean of the standard deviation of the TPM of each MAG for each plot⁴⁵. Spatial synchrony and population variability were calculated using the tempo and codyn R packages, respectively.

The contribution of stochastic and deterministic processes underlying the assembly of the biofilm community was quantified through a null-model analysis using the phylogenetic normalized stochasticity ratio (pNST). pNST was based on the β-nearest taxon index (βNTI)⁷². Briefly, the analysis compares the observed phylogenetic dissimilarity (β-nearest taxon index [βNTI])⁷² with that expected under null conditions⁷³. pNST metric ranges from 0, representing absolute determinism, to 100%, representing complete stochasticity. Statistical comparisons between treatments at the 2nd, 3rd, and 4th sampling dates were performed using a bootstrapping approach, with each dataset resampled 1000 times. P-values were adjusted for multiple comparisons using the false discovery rate (FDR) procedure, based on the number of comparisons at each sampling date. The analyses were performed with the pNST and nst.boot functions within the NST R package²⁴.

Since one of our primary goals was to assess the functional shifts underlying the response of biofilm to warming regimes, the calculation of functional entities was based on stress tolerance traits. Functional group richness (FR) was calculated as the number of functional entities (FE), which represent the number of MAGs sharing the same combination of stress tolerance traits. Because functional redundancy calculated as the ratio between the number of MAGs and FE may be biased, we computed functional over-redundancy (FOR) as proposed by Mouillot et al.²⁵. FOR represents the tendency of most MAGs to clump into a few highly abundant functional entities. FOR ranges from 0 to 1; low values of FOR indicate that most FEs contain the same number of MAGs, and high values indicate that most MAGs are contained in the richest FE. Finally, we computed functional vulnerability, which is the proportion of FEs represented by a single MAG. FEs and functional metrics were computed using the sp.to.fe and alpha.fd.fe functions, respectively, in the R package mFD⁷⁴. We used the number of stress-related gene modules in each MAG as a proxy for stress tolerance.

Statistical analysis

We used Linear Mixed Effects Models (LMEMs) to assess the effect of different temporal regimes of warming on taxonomic and functional diversity⁷⁵. The fixed part of the model included the “Temporal Regime of Warming” (5 levels: controls, fixed warming, fluct-s1+ fluct-s2+ fluct-s3) and “Time” (number days). To test the effect of fluctuating warming, “Temporal Regime of Warming” was partitioned into two planned contrasts: Control versus Fluctuating (Control vs Fluct. [fluct-s1+ fluct-s2+ fluct-s3] and Fixed versus Fluctuating (Fixed vs Fluct. [fluct-s1+ fluct-s2+ fluct-s3]). Each term and contrast were crossed with “Time”. The random part of the model included the variance in intercepts and slopes among plots \(({Time||plot})\). Since the main aim of the study was to assess the effects of varying warming regimes, the analysis involved the 2nd, 3rd, and 4th sampling date, for which treated plots received the same amount of warming. The analysis was centered on the 2nd sampling date so that deviations from the intercept (control) reflected the effects of the history of warming. In addition, to assess how warming history affected the recovery of biofilm following extreme temperatures, we performed contrasts equivalent to those of the main analysis, centered on the last sampling day. LMEM and contrasts were conducted using the glmmTMB and emmeans functions from the “glmmTMB” and “emmeans” R packages, respectively. P-values for contrasts at the last sampling date were adjusted using the false discovery rate (FDR) method. Model assumptions were assessed visually using plots of residuals vs fitted values, box plots of residuals vs. experimental conditions, and QQ plots of standardized residuals vs Normal quantiles generated with the “performance” R package⁷⁶.

A differential abundance analysis of MAG counts in fluctuating and fixed warming treatments relative to controls following temperature extremes was performed with the DESeq2 R package⁷⁷. To assess the trade-off between stress tolerance and growth rate, we examined the relationship between the number of gene modules related to stress-tolerance traits in each MAG and the minimal doubling time. Since the number of gene modules may also be influenced by genome size and MAG completeness, we included these variables as covariates in the model. Additionally, to account for variance in MDT, which changes with genome size, we applied a dispersion model (~log(genome_size)).

We also conducted an analysis to assess differences between control and extreme-only plots during the first phase of the experiment, before the application of extreme temperatures. The fixed effects included Treatment, Time, and their interaction, while the random effect was limited to plot ID (1|plot_id) to avoid convergence issues that arose when including Time as random slope.

Stability components

For each treatment, we computed four components of stability (sensitivity, resilience, temporal stability, and recovery) separately for biomass, dark and light yield using LMEMs^28,78. Treatment and time (number of days from the application of temperature extremes) were included in the fixed part of the model. Plot identity was included the random part. Time was also used as a covariate in the random part of the model to account for the lack of temporal independence of repeated observations within plots \(({{\rm{Time}}}||{{\rm{plot}}})\). Stability components were assessed over the five sampling dates following the application of temperature extremes pulses (from the 6th to the last sampling date). Controls were set as the baseline (Intercept) such that main effects (deviations between the Control and the treatments) reflected sensitivity after the imposition of extreme temperatures. Linear deviations between treatments and Control (treatment × time interactions) were used to quantify resilience. The following linear model was used:

Where y is the response variable (biofilm biomass) in the k-th plot and in the i-th treatment, \({\alpha }_{0}\) and \({\beta }_{1}\) are the intercept – centered at the first date after the imposition of extreme events – and the slope of controls; \({\alpha }_{i}\) is the main effect (sensitivity) of the i-th treatment, \({\beta }_{1i}\) is the deviation from the linear trend of the i-th treatment (resilience), \({a}_{0j}\) and \({b}_{1j}\) are the random effect intercept and slope of the j-th plot and \({\varepsilon }_{i,j}\) is the residual error.

Recovery was then estimated as the deviation between controls and each treatment at the last sampling date, using predicted values generated from LMEM (Eq. 1). Finally, temporal stability was quantified as the inverse of the standard deviation of the residuals. To derive measures of uncertainty for effects estimated through LMEM, we computed 1000 parametric bootstraps of each fitted model. Errors bars of sensitivity, resilience, temporal stability, and recovery are presented as the 95% confidence intervals (CIs).

To more accurately determine how a history of warming influences biofilm responses to extreme temperature events, we conducted a Before-After Control Impact (BACI) analysis. This analysis focused on two sampling dates before and after the temperature extremes. In our model, the fixed component encompassed both the treatment and the Before-After condition, as well as their interaction. The random component included the plot ID \((1|{{\rm{plot}}})\) and the Before-After condition, which was nested within the sampling date (Before-After:Time).

To control for potential artifacts, we conducted a separate analysis where the treatment was the sole covariate in the fixed part of the model. Here, the random structure included only the plot ID \((1|{{\rm{plot}}})\).

Reporting summary

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