Temperature variability interacts with mean temperature to influence the predictability of microbial phenotypes.

Despite their relatively high thermal optima (Topt ), tropical taxa may be particularly vulnerable to a rising baseline and increased temperature variation because they live in relatively stable temperatures closer to their Topt . We examined how microbial eukaryotes with differing thermal histories responded to temperature fluctuations of different amplitudes (0 control, ±2, ±4°C) around mean temperatures below or above their Topt . Cosmopolitan dinoflagellates were selected based on their distinct thermal traits and included two species of the same genus (tropical and temperate Coolia spp.), and two strains of the same species maintained at different temperatures for >500 generations (tropical Amphidinium massartii control temperature and high temperature, CT and HT, respectively). There was a universal decline in population growth rate under temperature fluctuations, but strains with narrower thermal niche breadth (temperate Coolia and HT) showed ~10% greater reduction in growth. At suboptimal mean temperatures, cells in the cool phase of the fluctuation stopped dividing, fixed less carbon (C) and had enlarged cell volumes that scaled positively with elemental C, N, and P and C:Chlorophyll-a. However, at a supra-optimal mean temperature, fixed C was directed away from cell division and novel trait combinations developed, leading to greater phenotypic diversity. At the molecular level, heat-shock proteins, and chaperones, in addition to transcripts involving genome rearrangements, were upregulated in CT and HT during the warm phase of the supra-optimal fluctuation (30 ± 4°C), a stress response indicating protection. In contrast, the tropical Coolia species upregulated major energy pathways in the warm phase of its supra-optimal fluctuation (25 ± 4°C), indicating a broadscale shift in metabolism. Our results demonstrate divergent effects between taxa and that temporal variability in environmental conditions interacts with changes in the thermal mean to mediate microbial responses to global change, with implications for biogeochemical cycling.

INTRODUCTION

The ocean has absorbed 50% of the carbon dioxide (CO2) emitted through anthropogenic activity in the Earth's biosphere and >90% of the heat trapped at the Earth's surface (Domingues et al., 2008). As a result, the global ocean has warmed at an average rate of ~0.062 ± 0.013°C per decade from 1900 to 2019 (Garcia‐Soto et al., 2021), and depending on the CO2 emissions scenario, it is anticipated surface waters will warm by at least twice that rate by 2100 under a low warming scenario (SSP1‐2.6) or by 4–8 times under a high warming scenario (SSP5‐8.5) (Arias et al., 2021). In addition to long‐term warming, the ocean has also experienced more frequent and intense marine heatwaves (anomalous warming events) over the last century (Oliver et al., 2018), with this trend projected to increase further under climate change (Frölicher et al., 2018).

The biological implications of a shifting baseline as well as increasing environmental variation are unknown (Kroeker et al., 2020). Our present understanding of the future distribution of species is based on climate‐envelope or ecological niche models, which forecast the availability of suitable environments for different taxa that are defined by their traits (Araújo et al., 2005; Dutkiewicz et al., 2013). The values of these traits are developed from field‐based observations and meta‐analyses that are largely derived from cultivation experiments under constant conditions (Dutkiewicz et al., 2015; Eppley, 1972). Although we are beginning to quantify the impacts of environmental variation compared with a change in mean condition (Cabrerizo et al., 2021; Gill et al., 2022; Pansch & Hiebenthal, 2019), such trait models have largely ignored the biological processes that influence species persistence under more variable conditions, and currently do not allow for organisms to dynamically respond to their environment.

Currently, it is unclear how selection for higher thermal optima (T opt) in response to ocean warming (i.e., increases in mean temperature) will affect species' capacities to adapt to changes in thermal variation. Plasticity, the phenotypic variation arising from the interaction of genotypes with their environment, is widely acknowledged as an important response to environmental change (Schlichting & Pigliucci, 1998). Theory predicts that plasticity should evolve in highly predictable environments, whereas reduced plasticity is expected in environments that fluctuate less predictably, because this lowers the risk that physiological responses do not match future selective pressures (Botero et al., 2015; Lande, 2009; Leung et al., 2020). There are a range of molecular, morphological, and physiological traits that are responsive to temperature, including body size (Daufresne et al., 2009), intracellular composition (Woods et al., 2003), the ratio of RNA to protein (Toseland et al., 2013), the use of heat‐shock proteins (HSPs) (Feder & Hofmann, 1999), and for photosynthetic microbes, photochemistry, carbon fixation and respiration, as well as changes in photosynthetic pigments and proteins like Rubisco (Baker et al., 2016; Schaum et al., 2018; Valentin & Mock, 2004). Temperature variability in marine environments shows a broad‐scale latitudinal pattern, with minimal seasonality in tropical environments and maximal in the temperate mid‐latitudes (Doblin & van Sebille, 2016; Locarnini et al., 2018). It follows that tropical species with highest T opt for growth and physiological functions have reduced tolerance to temperature variation due to their evolution under relatively stable temperatures (Qu et al., 2019; Thomas et al., 2012). Furthermore, the shape of thermal performance curves (TPC) means that fluctuations on the ascending portion of the TPC would lead to an acceleration in growth (Response A; see Figure 2), whereas fluctuations on the descending portion would cause a deceleration in growth (Response D) (Bernhardt et al., 2018).

FIGURE 2

Impact of temperature fluctuations on microbial eukaryotes. Carbon biomass normalised GPP (mean ± SD) and instantaneous growth rates (mean) of dinoflagellate strains in the downwards and upwards phase of fluctuating temperature treatments compared with growth in the stable control (n = 3). Response A is characterised by growth in the downward/upward phase being lower/higher than control (mean) temperature. Response D is characterised by growth in both the downward and upward phase being lower than the stable control. Note different x‐axis scales between plots.

Here, we explicitly test how temperature variability interacts with mean temperature to influence growth and phenotypic plasticity among cosmopolitan tropical and temperate eukaryotic microbes. We used two species of the same genus, and two strains of the same species, each with a different thermal history, to investigate how universal their responses were. Strains included the potentially harmful dinoflagellate Coolia palmyrensis and Coolia malayensis with similar thermal traits (T opt, niche breadth) even though they were isolated from tropical and temperate latitudes, respectively (Larsson et al., 2019), and two strains of the tropical dinoflagellate Amphidinium massartii with different thermal traits; one that was maintained at the standard growth temperature of 25°C (referred to as the control strain [CT]) and another that had undergone directional high‐temperature selection at +5°C for 3.25 years (>500 generations) (referred to as the high‐temperature strain [HT]) (Baker et al., 2018). Temperature treatments were designed to fluctuate with a period of 2 days at different amplitudes (0 control, ±2, ±4°C) around different means (20, 25, 30°C) that were below or above the estimated T opt of each strain (Figure 1; Table S1). Exposure to regular temperature fluctuations led to a universal reduction in average population growth rate, but there was a greater cost to growth for temperate C. malayensis and A. massartii HT (strains with narrower thermal niche breadth). At suboptimal mean temperature, traits were linearly correlated to cell volume. However, at supra‐optimal mean temperature, novel trait combinations developed, leading to greater phenotypic heterogeneity. Furthermore, growth declined as expected in the warm phase but also declined during the cool phase of fluctuations, indicating that the TPC does not reliably predict performance at supra‐optimal mean temperature. Together, these empirical data show that phenotypic traits become less predictable at mean temperatures above T opt and thermal specialisation (narrower niche breadth) may come at the expense of lowered tolerance to temperature variability.

FIGURE 1

Design of temperature treatments. Temperatures fluctuated regularly every 2 days around a mean throughout an initial acclimation period involving six growth cycles (20–28 days). Cells were then harvested under the final downward (cool) and upward (warm) phase to quantify phenotypic traits.

MATERIALS AND METHODS

Microbial eukaryote taxa and temperature selection environments

Experiments were performed using dinoflagellates isolated from tropical and temperate latitudes (Table S1). Tropical A. massartii (CS‐259) was provided by S. Murray after it was isolated from Kurrimine Beach and described in Murray et al. (2012). It was grown under temperature selection for 39 months at two different temperatures (25 and 30°C; CT and HT, respectively), after which the HT strain had become high‐temperature specialized, with a +2.67°C shift in the T opt and a +14.92°C increase in the T min, with the overall effect of narrowing the thermal niche (Baker et al., 2018). A second tropical dinoflagellate C. palmyrensis and a temperate dinoflagellate Coolia malayensis isolated within months of each other in 2014 (Larsson et al., 2019) were also used, allowing us to examine whether thermal history influences responses of strains to thermally variable environments.

A. massartii was cultivated in coastal seawater medium with modified f/2 enrichment, lacking silicic acid (Guillard & Ryther, 1962). Coolia spp. was cultivated in coastal seawater medium with modified K nutrients (Litaker et al., 2009). All cultures were maintained under a light intensity of 100 μmol photons m−2 s−1 with a light–dark cycle of 12:12 light: dark (LD) in either 25°C (CT, Coolia spp.) or 30°C (HT).

Experimental setup

Thermal regimes were designed to test the effect of regular thermal fluctuations in comparison with stable control treatments with the same mean temperature. Control treatments had thermal variability of <0.1°C over the experimental duration. The amplitude of temperature fluctuations was ±2 or ±4°C with temperatures in the variable treatments cycling every 2 days between temperature endpoints (Figure 1). Similar frequencies of short‐term temperature variations are observed in coastal waters worldwide (Aguirre et al., 2021; Dai et al., 2009; Leinweber et al., 2009; Shaw et al., 2012). Such temperature variations are especially pronounced in coastal regimes due to the effects of atmospheric weather fluctuations on a shallow water column along with vigorous advection, mixing and upwelling (Aguirre et al., 2021; Leinweber et al., 2009). Mean temperatures for experiments were determined based on two considerations: (1) seasonal sea surface temperatures in tropical eastern Australian waters (AIMS, 2009): 20°C (winter; July) and 30°C (summer; January) and (2) the long‐term culturing temperature of isolates (Table S1). For CT ad HT strains, this equated to reciprocal temperature transplants at 25 and 30°C.

To determine the effect of amplitude in variable thermal regimes, all taxa were exposed to 25 ± 2 and ±4°C, relative to the stable control—that is, temperature varied from 23°C to 27°C or from 21 to 29°C (Table S2). To examine the effect of thermal variation centered around different mean temperatures, A. massartii CT and HT strains were exposed to 30 ± 4°C (26 to 34°C), and C. palmyrensis and C. malayensis were exposed to 20 ± 4°C (16 to 24°C; Table S2). Given the TPC of Coolia spp., fluctuations around 20°C were below T opt and fluctuations around 25°C were supra‐optimal; for Amphidinium, fluctuations around 25°C were around T opt, and those around 30°C were supra‐optimal.

Experiments were conducted in a temperature‐controlled room, where a system of water baths was used to maintain incubation temperatures. Temperatures in the variable treatments were adjusted every 2 days using thermostatically controlled cooling and heating systems (Julabo GmbH), with temperature shifts occurring over 60 min. Experimental microcosms consisted of 500 ml capped polycarbonate flasks filled with 400 ml of culture, illuminated from above with LED light banks (Schenzen Cidly Group) at 150 μmol photons m−2 s−1 in a 12:12 light: dark cycle.

Semi‐continuous dilution culturing methods were used to acclimate cells to variable temperature treatments (Figure S2) and avoid excessive biomass accumulation and consequent nutrient limitation or self‐shading. At the end of every dilution cycle, each bottle was diluted individually back to the initial cell concentration for that bottle. Thus, in this semi‐continuous “turbidostat” method, the growth rate in each treatment determined the dilution rate. After 6 cycles (~32 days), cells were subsampled for physiological trait quantification at the end of the 2‐day cool exposure period (downward phase) and then 2 days later, at the end of the warm exposure period (upward phase) (Figure 1). All experiments used three replicate bottles in each treatment.

Trait quantification

Biomass, relative cell size estimates, and specific growth rates

Biomass was estimated every second day using in vivo chlorophyll‐a (Chl‐a) fluorescence to determine semi‐continuous dilution volumes based on calculated specific growth rates. These in vivo Chl‐a estimates were later validated with microscopic or flow cytometric cell counts; growth rates calculated with these two methods were found to be closely comparable.

Cell abundance of A. massartii was determined with flow cytometry. One mL of culture was first fixed in paraformaldehyde solution (1% v/v final concentration), snap frozen in liquid nitrogen and stored at −80°C. Prior to enumeration, samples were thawed quickly in a 35°C water bath, vortexed and analysed with an Influx flow cytometer (BD Biosciences) under blue (488 nm laser) excitation. Yellow‐green fluorescent beads were added in known concentrations for volumetric and fluorometric internal calibration. Events were triggered on FSC for enhanced sensitivity; populations were isolated on cytograms of FSC and chlorophyll‐a fluorescence (692/40 nm). Cell count samples for the Coolia experiments were preserved in Lugols (1% final solution) and were then enumerated under the microscope using a Sedgewick Rafter Chamber.

Cell size

A. massartii cell size was estimated using flow cytometric forward scatter, relative to bead standards (Marie et al., 1999). For Coolia spp., Lugols‐preserved subsamples were aliquoted into 24‐well glass bottom plates, and cells were measured with an automated imaging system (IN CELL 2200, GE Healthcare) using a 20× objective. Images were analysed with proprietary software to obtain major axis length, minor axis length and form factor (i.e., cell roundness, varying between 0 and 1 with 1 being a perfect circle). Biovolume was calculated based on an ellipsoid shape using minor axis as width and major axis as length, according to Hillebrand et al. (1999).

Chlorophyll and elemental composition

Chlorophyll‐a samples were extracted using 90% acetone for 24 h in a −20°C freezer in the dark and measured on a calibrated fluorometer (Turner designs Model 10 AU) (Fu et al., 2005). Particulate organic carbon, nitrogen and phosphorus (POC, PON, and POP, respectively) samples were filtered onto pre‐combusted GF/F filters and dried at 60°C, then stored at room temperature. PON and POC was measured using an elemental combustion system (Model 4010, Costech) (Hutchins et al., 1998). POP samples were combusted at high temperature (500°C) to convert the organic P into orthophosphate, and then measurements were performed using the standard colorimetric molybdate method (Strickland & Parsons, 1968).

Primary productivity (14C uptake)

To estimate primary productivity, carbon fixation rates were measured using 14C‐labeled bicarbonate in small volume incubations as described in Doblin et al. (2011). Specifically, radiolabelled NaH14CO3 (stock solution 1.85 × 107 Bq) was added to 5 ml of culture in clear glass tubes (1.5 μCi per tube) and incubated at the treatment temperature for 60–80 min under the growth irradiance. Activity in the samples was determined by removing a 100 μL aliquot and placing it into 5 ml of refrigerated 0.1 M NaOH, adding 10 ml scintillation fluid (Ultima Gold™, PerkinElmer) and shaking before counting using a liquid scintillation counter (Packard TriCarb 2900TR). Following incubation, tube contents were acidified with 250 μl 6 M HCl and shaken on an orbital shaker for 12 h to remove unfixed 14C. Scintillation fluid (10 ml Ultima Gold, Perkin Elmer) was then added to each sample, vigorously shaken and left for 1 h before counting. Counting time was set to 5 min so that counts were within a 5% counting error.

RNA sampling, extraction, and analysis

Cells were sampled between hour 1 and 2 of the light period to limit diel effects, and immediately filtered onto 5 μm 25‐mm polycarbonate membranes before being flash frozen in LN2 and stored at −80°C. Due to resource limitations, only control and ±4°C fluctuation samples from CT, HT, and C. palmyrensis were analysed further (Table S2). RNA extraction was performed based on methods described in Verma et al. (2019). Briefly, preheated Trireagent (Ambion) and bead beating was used for cell lysis. RNA was purified using RNeasy Mini Kit (Qiagen) according to manufacturer's instructions. Any residual DNA was removed via the TURBO DNA‐free™ Kit (Life Technologies) according to manufacturer's instructions and RNA was then stored at −80°C. The RNA purity, quantity, and integrity were assessed using a Nanodrop ND‐1000 (Thermo Scientific) and 2100 Bioanalyser Nano and Pico chips (Agilent Technologies) depending on RNA yields.

Following RNA extractions, sequencing libraries were prepared for 2–3 replicates per sample using the Illumina TruSeq stranded mRNA sample prep following the manufacturer's instructions in two batches. The first batch comprised all libraries from Amphidinium 30 ± 4°C samples. Those samples yielded relatively low RNA quantities (on average 491 ng), hence the libraries were prepared using a higher number of PCR cycles (15 cycles). The second batch comprised the libraries from all remaining samples, prepared using 12 PCR cycles. All libraries (71) were sequenced together using one Illumina NovaSeq S1 flow cell (paired‐end, 2 × 100 bp reads) yielding 23 to 58 million read pairs per sample (30 million on average). Raw sequencing data has been deposited to the NCBI Sequence Read Archive under Bioproject PRJNA819215.

Quality trimming, assembly, clustering, differential expression analysis

Sequencing reads were quality trimmed using trimmomatics (Bolger et al., 2014), using the default settings (ILLUMINACLIP:$TRIMMOMATIC_DIR/adapters/TruSeq3‐PE.fa:2:30:10 SLIDINGWINDOW:4:5 LEADING:5 TRAILING:5 MINLEN:25). Quality trimmed reads were assembled into transcripts using Trinity version 2.6.6 (Grabherr et al., 2011; Haas et al., 2013). For C. palmyrensis, reads from all samples and replicates were used to generate one combined assembly. Similarly, for A. massartii, because both strains originate from the same clonal culture, reads from all samples and replicates from the CT and HT strains were used to generate one combined Amphidinium assembly. In each assembly, similar transcripts were clustered and only the longest transcript retained using CD‐HIT‐EST (word size of 8 and sequence identity threshold of 0.9) (Weizhong & Adam, 2006). Transcriptome completeness was estimated using BUSCO (version 3.0.2) and the eukaryote database of conserved single copy genes (eykaryota_odb9, creation date: 2016‐11‐02) (Simão et al., 2015).

Transcript abundance was calculated using kallisto (Bray et al., 2016) within the Trinity package. Differential expression analysis was done using edgeR (Robinson et al., 2010) within the Trinity package and was performed separately for each species and experiment, using a false discovery rate (FDR) of <0.05 and a minimum twofold change in abundance. To identify differentially expressed transcripts between the variable and the stable control treatments, differential expression analysis was performed on each sample of the variable treatments against both samples from the respective control treatment; only transcripts that were significantly DE in comparison with both control samples were considered DE between variable and control treatment.

Functional annotation of transcripts, including Gene Ontology (GO) annotation, was done using Blast2GO version 5.2.5 (Gotz et al., 2008) based on BLASTX searches against the NCBI non‐redundant protein database (nr) (standalone BLAST+ version 2.7.1 [Camacho et al., 2009]) and InterProScan searches (version 5.30‐69 [Jones et al., 2014]). Annotation of photosynthesis and stress‐related transcripts was manually confirmed. Enrichment testing of biological process (BP) GO categories within sets of differentially expressed transcripts was done within Blast2GO using Fisher's exact test with a FDR of 0.05, corrected for multiple hypothesis testing (Benjamini & Yosef, 1995). Enriched BP GO categories were further reduced to only the most specific terms (more general GO terms on higher levels in the Directed Acyclic Graph [DAG] were removed).

Statistical treatment of data

Statistical analyses were performed using IBM spss Statistics 24. Prior to analysis data was checked for normality and homoscedasticity. Data were analyzed by ANOVA using strain, mean temperature, and amplitude as fixed factors. Differences in multi‐variate phenotypes were tested using ANOSIM (PRIMER‐E v.6). Draftsman plots were used to visualise correlations among traits before the similarity matrix using Bray‐Curtis distance was constructed.

RESULTS

Effects on growth

Following an acclimation period (Figure 1 and Figure S2), average population growth rate under fluctuating conditions was significantly lower than under stable conditions in all four dinoflagellate strains (Table 1; ANOVA Type III, F [1,23] 23.93, p < .001), irrespective of mean temperature. Thermal variability caused the greatest negative impact (>70% decline) to both the HT and CT strains of A. massartii when the temperature fluctuated ±4°C around the supra‐optimal (30°C) compared with the standard mean (25°C) (ANOVA Type III, F (2,22) = 132.561, p < .001). The amplitude (±2 vs. ±4°C) of the temperature fluctuation (tested at 25°C for all taxa; see Section 2 and Table S2) was significant for A. massartii HT (ANOVA type III, F [2,6] = 54.349, p < .001) and tropical C. palmyrensis (ANOVA Type III, F [2,6] = 40.776, p < .001; Table 1), both showing a greater decline in the average population growth with increasing amplitude. However, the decline in growth was similar in the ±2 and ±4°C fluctuations for A. massartii CT and the temperate C. malayensis population, which may be due to larger between replicate differences.

TABLE 1

Growth rates (mean ± SD; day−1) for four dinoflagellate strains grown at different mean temperatures with different amplitude of regular (2‐day) temperature fluctuations (n = 3)

	20°C mean	20°C mean	25°C mean	25°C mean	25°C mean	30°C mean	30°C mean
	Stable control	±4°C	Stable control	±2°C	±4°C	Stable control	±4°C
Amphidinium massartii CT			0.630 ± 0.019a	0.540 ± 0.010ab	0.513 ± 0.059b	0.686 ± 0.016	0.199 ± 0.012*
A. massartii HT			0.527 ± 0.027a	0.420 ± 0.009b	0.364 ± 0.018c	0.557 ± 0.027	0.107 ± 0.015*
Coolia palmyrensis tropical	0.264 ± 0.010	0.197 ± 0.007*	0.412 ± 0.008a	0.370 ± 0.005b	0.346 ± 0.012c
C. malayensis temperate	0.245 ± 0.002	0.163 ± 0.017*	0.359 ± 0.012a	0.269 ± 0.008b	0.295 ± 0.032b

Note: Asterisks show that growth rates under fluctuating temperature are significantly lower than those at stable temperature (ANOVA, p < .05), with letter superscripts indicating differences among different amplitudes at a mean temperature of 25°C.

Growth of all strains responded dynamically to within treatment temperature changes, with instantaneous growth rates diverging at temperature fluctuation endpoints (i.e., cool/downward and warm/upward phases; Figure 2). Reduced growth during downward fluctuations and increased growth during upward fluctuations (with intermediate growth under stable mean temperature) was anticipated at mean temperatures below T opt in the ascending portion of the TPC (Bernhardt et al., 2018). This response was termed Response A (Figure 2). In contrast, Response D occurred at mean temperatures above T opt in the descending portion of the TPC and was characterised by reductions in growth rate under both upward and downward temperature fluctuations relative to the stable temperature treatment, indicating that cells were redirecting resources away from cell division (Figure 2).

The average population growth rate of Coolia species under fluctuating temperatures was less impacted than Amphidinium (Figure 2; Table 1), likely due to longer generational times (2.6 days in Coolia spp. vs. 1.7 days in Amphidinium: Table 1) that exceeded the 2‐day frequency of fluctuations. At the mean temperature of 20°C (3.6 to 4.5°C below T opt; Table S1), instantaneous growth rates of Coolia spp. under ±4°C fluctuations diminished by −96.0% ± 10.5% (mean ± SD; tropical C. palmyrensis) and −110.4% ± 1.2% (temperate C. malayensis) relative to stable controls in the downwards phase of the fluctuation (16°C) and increased to 64.8% ± 18.3% and 43.6% ± 13.2% in the upwards phase (24°C)—that is, the growth dynamics of both species were characteristic of Response A (Figure 2). However, under ±4°C temperature fluctuations at the mean temperature of 25°C (0.5–1.4°C above T opt; Table S1), the tropical C. palmyrensis strain exhibited Response A (instantaneous growth rates declined by −39.4% ± 5.6% at 21°C and increased by +8.1% ± 1.2% at 29°C) whereas the temperate C. malayensis showed Response D, characterised by a reduction in instantaneous growth in both the downwards (−36.3% ± 0.2%) and upwards phases (−9.4% ± 13.3%) of the 25 ± 4°C fluctuation relative to stable controls (Figure 2).

Both strains of A. massartii exhibited similar growth patterns under thermal fluctuations. At the mean temperature of 25°C with ±4°C fluctuation, CT and HT exhibited Response A. Growth of the CT strain (T opt 23.94°C) declined by −58.4% ± 12.7% in the downwards phase (21°C) and increased by 24.1% ± 19.9% in the upwards phase (29°C) relative to stable controls (Figure 2). The HT strain (T opt 26.61°C) showed an even greater decline in growth at 21°C in the downwards phase (−88.4% ± 4.1%) compared with the CT strain (Figure 2), and a similar increase in growth (+27.0% ± 2.7%) at 29°C in the upwards phase. At mean temperature of 30°C (+3.39 to 6.06°C above T opt for both strains; Table S1), both CT and HT exhibited Response D but the decline in growth was 10% greater in the HT compared with the CT strain (−76.7% ± 5.5% vs. −66.4% ± 3.9% at 26°C and −83.7% ± 3.2% vs. −75.6% ± 1.7% at 34°C, respectively; ANOVA F [1,9] = 14.508, p = .004), resulting in a significantly lower overall HT growth rate (Figure 2; Table 1).

Effects on phenotypic traits

There were significant trait value changes between the downward and upward phases of the fluctuations, indicating phenotypic plasticity of all strains (Data Table). As part of Response A, cell volume was directly proportional to growth (Figure 3; Table S3), most clearly illustrated by the tropical C. palmyrensis in the 20 ± 4°C treatment. When temperatures approached 16°C, growth slowed (i.e., cell division ceased), and cell volume was maximal, as was elemental carbon (C), nitrogen (N), and phosphorus (P) content. Conversely, as temperatures approached 24°C, growth increased, and cell volume diminished alongside cellular C, N, and P. Cellular chlorophyll‐a (Chl‐a) content and gross primary productivity (GPP; carbon fixation measured over 1–2 h) showed poor correspondence to cell volume, but C:Chl‐a was positively correlated. The same pattern was observed for both Amphidinium strains when temperatures fluctuated around 25°C (Figure 3).

FIGURE 3

Cell volume–dependent growth. Temperate and tropical Coolia spp. strains (green and orange symbols, respectively) at 20°C (a), all strains at 25°C (b) and Amphidinium massartii control temperature and high‐temperature strains (blue and red symbols, respectively) at 30°C (c). Lines reflect significant linear regression (p < .05). Size of symbol reflects magnitude of fluctuation amplitude, whereby largest symbols are ±4°C treatments, smallest symbols are stable treatments and intermediate symbols are ±2°C treatments. Upward and downward facing symbols depict warm and cool phases of temperature fluctuation treatments, respectively. Each symbol represents a distinct biological replicate.

Response D however, characterised by relatively low growth rates at the downward and upward phases of temperature fluctuations around a supra‐optimal mean, led to new trait combinations that were most apparent in Amphidinium (Figure 4). In the 30 ± 4°C regime, A. massartii cells became significantly larger but a distinct elemental composition was evident for the HT strain which had higher N per biomass and hence lower C:N compared with the CT population (Figure 4). In the 25 ± 4°C treatment close to T opt, Amphidinium cells were smallest, most pigmented and had maximal rates of carbon fixation (GPP; gross primary production). Amphidinium phenotypes diverged between stable and fluctuating treatments (ANOSIM Global R = .330, p = .001), at specific temperatures (ANOSIM Global R = .788, p = .001) except for 25 and 29°C (pairwise tests) as well as between CT and HT strains (ANOSIM Global R = .117, p = .019). Coolia phenotypes appeared less plastic to the imposed thermal treatments (Figure 4), with tropical and temperate species responding similarly (ANOSIM, R = .056 p = .113) with phenotypes only distinct between different temperatures (Figure S5; ANOSIM Global R = .678, p = .001).

FIGURE 4

Viable multitrait phenotypes of Coolia spp. and Amphidinium massartii. Resultant trait combinations in control and fluctuating temperature centred around a mean of 20°C (Coolia) and 30°C (Amphidinium). nMDS input variables include growth, cell size, cellular carbon, nitrogen, phosphorus content, C:N, N:P, chlorophyll‐a content, and gross primary productivity. Vectors show the traits driving differences between phenotypes (Pearson correlation >0.7). Symbols represent tropical Coolia palmyrensis (orange) temperate Coolia malayensis (green), A. massartii CT (blue) and HT (red) strains. Symbol directions (upward or downward) and size are the same as Figure 3.

Effects on metabolism

To help uncover how changes in the temperature environment impacted different strains, we examined carbon fixation (i.e., GPP) in relation to growth rates as well as gene expression at downward and upward phases relative to stable mean temperature controls. In Coolia spp., GPP (normalised to carbon biomass) was proportional to growth at 20 ± 4°C, increasing with growth rate in the upward phase, and decreasing with growth rate in the downward phase of fluctuations (Figure 2). However, when temperature fluctuated at 25 ± 4°C, GPP was a factor of ~2 higher and remained relatively constant in both downward and upward phases, with instantaneous growth rates that were similar to the 20 ± 4°C treatment (Figure 2). The pattern was similar in Amphidinium, but because growth did not increase in the downward phase at supra‐optimal mean temperature (Figure 2; Table 1) we interpret this to mean that fixed carbon was being directed away from cell division. Notably, biomass normalised carbon fixation rates were lower in the HT versus CT strain except in stable temperature controls (Figure 2).

In Coolia spp., growth rates were lowest in the downward phase of fluctuations at mean temperatures of both 20 and 25°C for both the tropical and temperate strains (Figure 2). Even though the tropical C. palmyrensis strain grew extremely slowly at 16°C (downward phase of 20 ± 4°C) (0.011 ± 0.028 day−1; ~97% slower than its maximum), there was relatively little downregulation of gene expression (Figure S1). In comparison, a wide range of transcripts encoding photosynthetic functions were upregulated at both 21 and 29°C (25 ± 4°C), consistent with higher carbon fixation rates than in the 20 ± 4°C regime (Figure S3). When exposed to 29°C (upward phase of 25 ± 4°C), instantaneous growth rates of C. palmyrensis were faster than in the stable mean control temperature of 25°C (0.445 ± 0.005 day−1 vs. 0.412 ± 0.008 day−1, respectively) and were accompanied by large transcriptional changes (Figure S1). Although the number of downregulated transcripts was large, many of them could not be annotated and thus encoded for unknown functions. In contrast, 43 BP GO categories for various metabolic functions were enriched in the upregulated transcripts (Figure S4), including 28 transcripts involving photosynthesis functions (Figure S3B), indicative of a broad‐scale shift in metabolism.

In general, more transcripts were differentially expressed in the A. massartii HT strain than the CT strain in both the 25 ± 4°C and 30 ± 4°C regimes, suggesting it was undergoing more physiological adjustments (Figure S1). At 21°C (downward phase of 25 ± 4°C), the HT strain had a considerably slower growth rate at 0.061 ± 0.022 day−1 (~90% slower than its maximum) relative to the CT strain (0.262 ± 0.080 day−1), consistent with it having a higher T min (14.94 vs. 0.02°C for CT; Table S1). It also downregulated considerably more transcripts and functions (2008 transcripts within 27 enriched BP GO categories compared with 94 transcripts and 7 enriched BP GO categories in CT; Figure 5a, Figure S1), including 26 transcripts pertaining to the photosynthetic light harvesting complex (chlorophyll A–B binding proteins) (Figure 5b). DNA repair enzymes and HSPs and related chaperones were also encoded among the downregulated transcripts in the HT strain at 21°C (25–4°C) (Figure 5C), suggesting a broad‐scale shift in metabolism. In the upwards phase of each temperature regime, all strains were actively dividing, despite being exposed to absolute temperatures approaching their thermal maximum (T max) (Table S1). Under periodic exposure to 34°C (upwards phase of 30 + 4°C for Amphidinium only), both the HT and CT A. massartii strains showed a large number of differentially expressed transcripts relative to the 30°C stable mean control (2523 and 3374 for HT and CT, respectively), with the majority of those transcripts upregulated at 34°C (78% and 76% for HT and CT, respectively) (Figure S1). In both Amphidinium strains, transcripts encoding functions involved in DNA recombination and integration were enriched in the upregulated transcripts, indicative of genomic changes at 34°C (Figure S4).

FIGURE 5

Transcriptional plasticity of dinoflagellates. Number of biological process gene ontology categories enriched in differentially expressed transcripts (a), differentially expressed transcripts encoding photosynthesis (b) and stress functions (c) for A. massartii control temperature and high‐temperature strains growing at fluctuating temperatures (25°C ± 4 (21 and 29) and 30°C ± 4 (26 and 34)) relative to stable controls. “down/up” refers to downregulated and upregulated differentially expressed transcripts, respectively.

Together, these observations suggest that variable temperatures approaching upper thermal limits led to increased metabolism, despite contrasting effects on growth, and that in fluctuating environments, cells may tolerate periodic exposure to temperature extremes with limited signs of stress.

DISCUSSION

A greater understanding of the consequences of thermal variation are needed to make predictions about the current and future growth, abundance, and distributions of species in the natural environment. Climate change models not only forecast higher mean temperatures in the future ocean but also that the magnitude, frequency, and extremes of temperature variability will increase (Boyd et al., 2016; Frölicher et al., 2018; IPCC, 2014, 2018; Thornton et al., 2014; Vasseur et al., 2014). Here we show that regular 2‐day temperature fluctuations of the type often observed in coastal waters (AIMS, 2009; Dai et al., 2009; Leinweber et al., 2009; Shaw et al., 2012), can cause significant demographic and biogeochemical impacts on microbial eukaryotes, and that phenotypic trait values diverge at a similar absolute temperature depending on previous thermal history. Furthermore, the existence of Response D (decline in growth under both upwards and downwards shifts in temperature from a supra‐optimal mean) indicates that predicting responses is not as straight forward as understanding the shape of TPCs.

Our expectation was that average population growth rates would decrease under variable temperature treatments because cells spent less time within their T opt range, and that the decline would be proportional to the amplitude of temperature fluctuation because increasingly variable regimes impose an energetic cost (Putnam & Edmunds, 2011). Given that TPCs are non‐linear (Huey & Stevenson, 1979), with an accelerating phase of growth at temperatures below T opt and a decelerating phase of growth at temperatures greater than T opt, we anticipated Response A, where instantaneous growth rates decreased during the downward phase and increased during the upwards phase of fluctuations (Bernhardt et al., 2018). However, there was a net loss of fitness when temperature fluctuated around a supra‐optimal mean, as demonstrated by the non‐reversibility of growth (i.e., Response D) in temperate C. malayensis and A. massartii CT and HT—that is, no increase in growth rate during the downward phase as temperature approached T opt. This indicates there is an interaction between temperature variability and mean temperature that fundamentally alters cellular processes, with cascading effects on growth and phenotypic traits. One direction for future work would be to determine the threshold temperature above T opt and associated fluctuations where phenotypic responses are no longer reversible, and what other factors affect the temperature at which this thermal tipping point occurs.

Although there is well‐established theory that overall phenotypic plasticity increases when populations evolve in predictably variable environments, and some empirical evidence that environmental fluctuations with fixed period but random amplitude select for plastic individuals (Botero et al., 2015; Lande, 2009; Schaum et al., 2016), the combinations of traits that phytoplankton use to increase plasticity is less well known. The elemental content of C, N, and P was proportional to cell volume which flexed depending on growth rate at temperatures below T opt, demonstrating its utility as a “master trait” (Litchman & Klausmeier, 2008). Given that C, N, and P make up most microbial biomass (Finkel et al., 2016), changes in these traits with cell size were expected. However, pigment content and GPP were not related to cell size, and at temperature extremes it was evident that fixed carbon was being redirected into non‐growth components—for example, respiration or dissolved organic carbon release or both. Our observations also showed that the A. massartii HT strain diverged from its CT descendant, showing higher N content than at the stable mean temperature of 30°C, consistent with its increased N demand associated with its high‐temperature adaptation (Baker et al., 2018). Furthermore, C:N and N:P ratios were related to cell volume but showed contrasting patterns in Amphidinium versus Coolia. There was also a common pattern of increasing heterogeneity among phenotypes at temperatures exceeding T opt. Although many ocean biogeochemical models include correlations between traits (e.g., those that incorporate flexible C:N:P stoichiometry), many overlook how these correlations may change with adaptation to climate change. Collapsing multi‐trait phenotypes, derived from experiments such as this, into two dimensions using multi‐variate methods may provide a pathway for integrating plasticity and evolution into ocean biogeochemical models (Argyle et al., 2021; Walworth et al., 2021).

Another aspect of plasticity revealed by transcriptomic analyses in this study was that microbial eukaryotes can be in physiologically different states while exhibiting similar growth rates. The ability to down‐ and up‐regulate metabolism whilst maintaining a consistent growth rate was demonstrated by the tropical C. palmyrensis strain fluctuating between 21 and 29°C (25 ± 4°C; Figure 2 and Figure S4). Alternatively, broadscale downregulation may be a conservative strategy that protects cells against damage when approaching temperature limits (Kültz, 2005) but requires cells to rapidly respond to changes in their external conditions. The response of C. palmyrensis to the sub‐optimal temperature at 16°C (downward phase of 20 ± 4°C) involved cell persistence and cessation of cell division (Figure 2), but no major restructuring of metabolism (Figures S3 and S4). Importantly, our gene expression analyses showed that temperatures of upper and lower thermal growth limits determined under stable conditions (Table S1) were not necessarily stress‐inducing or fatal for cells that only periodically experienced these extreme temperatures. Temporary respite from high or low temperature evidently allowed populations to limit cumulative stress, as was shown by Schaum et al. (2018) for a marine diatom in fluctuating temperature regimes. These observations are in line with field data that show similar composition of dinoflagellates within phytoplankton communities in the tropical Pacific Ocean despite a clear environmental gradient in temperature and nutrient limitation (Cohen et al., 2021), suggesting that regulation of gene expression is critical to plasticity in situ. However, thermal tolerance within any single taxon could be surpassed during marine heatwaves (temperatures above a seasonally varying 90th percentile for at least 5 days at a particular location; Hobday et al., 2018). Indeed, in mixed natural phytoplankton communities from the San Pedro Ocean Time‐series station exposed to mean temperatures that exceeded multiyear upper thermal limits and fluctuated ±4°C, there was a complete restructuring of the community (Kling et al., 2020).

At the outset of this study, we hypothesised that the thermal histories (i.e., differing isolation locations or laboratory selection regimes) and higher T opt of A. massartii HT and tropical C. palmyrensis would lead to greater vulnerability to temperature variability. However, this was not the case, with narrower thermal niche breadth the common factor among strains that had the most detrimental impacts to growth under fluctuating temperature.

The mean temperatures used in our experiments were 20, 25, and 30°C, representative of temperate, sub‐tropical and tropical ocean waters, respectively, and therefore relevant to species growing in a large proportion of the global ocean (Locarnini et al., 2018). Fluctuation regimes were designed to have regular intervals (2 days) that would not be biologically anticipated through diel light cues (Ottesen et al., 2014; Vislova et al., 2019). In general, A. massartii populations were more responsive to temperature shifts, moderating growth rates up and down when temperatures fluctuated. Given that fluctuations occurred every 2 days, this exposed Amphidinium cells to temperatures that allowed for cell division, whereas for Coolia spp., temperature exposures were within a single generation. This suggests that different generational times could set the pace of acclimatization relative to the frequency of environmental change—that is, fast‐growing taxa with short generation times could potentially alter their cell physiology at similar time scales to environmental fluctuations. Among freshwater phytoplankton taxa Fey et al. (2021) found greater support for gradual thermal acclimation that lagged behind changes in temperature, rather than instantaneous adjustments. Furthermore, mechanisms such as nutrient storage could allow large cells to “ride out” changes in the external environment. Considering that growth rates are intrinsically related to cell size (Litchman et al., 2007), our results potentially provide a unifying framework with which to test the implications of environmental fluctuations on microbial taxa spanning orders of magnitude in size.

CONCLUSION

With the ocean projected to become warmer and more variable, investigations that connect aspects of environmental variation with organismal plasticity and evolution are urgently needed (Collins et al., 2020). Here we show clear evidence that microbial eukaryotes diverge in their responses to temperature fluctuations depending on whether mean temperature is above or below T opt, with the magnitude of impact being mediated by their degree of thermal specialisation. Our selection of strains allowed comparison of dinoflagellate taxa with similar thermal traits isolated from different latitudes (Coolia spp.) as well as a comparison of strains from the same ancestor that were maintained at different temperatures (A. massartii), both of which showed clear differences linked to narrower thermal niche breadth not T opt as we hypothesised. The implications are that species with narrower TPCs under a warming baseline could be more vulnerable to changes in temperature variation, with cascading impacts to global biogeochemical cycles and climate feedbacks.

FUNDING INFORMATION

This research was funded by the Australian Government through the Australian Research Council Discovery Scheme (project no. DP140101340 awarded to MD, SM, and DH and DP180100054 awarded to MD). FF and DH were funded by a University of Technology Sydney Distinguished Visiting Professor Scheme. Additional support was provided by U.S. National Science Foundation grant OCE 1638804 to DH and OCE 1851222 to FF and DH. KB, ML were each supported by an Australian Government Training Program Scholarship.

CONFLICT OF INTEREST

MD is a member of the board of directors for the Sydney Institute of Marine Science.

Appendix S1

Click here for additional data file.