Introduction

Phytoplankton, which lie at the base of most marine food webs, play crucial roles in global biogeochemical cycles and contribute 50% of the primary production of the biosphere [1, 2]. These organisms face multiple environmental changes, including rising sea surface temperature (ocean warming) and increasing CO2 levels in the context of global change [3]. It is predicted that phytoplankton productivity will decrease, and the community structure will become increasingly unstable, in the future high CO2 and warming ocean [4, 5]. However, such models ignore the evolutionary responses of phytoplankton to the environmental changes of high CO2 and/or warming [6, 7]. Such predictions in the absence of evolutionary responses would add to the uncertainty of the projected changes in phytoplankton productivity, community structure, and diversity. On the one hand, while it has been demonstrated that marine phytoplankton can promote swift evolutionary responses to adapt to high CO2 or warming alone [8,9,10], little is known regarding the molecular mechanisms underpinning the evolutionary responses of phytoplankton to a combination of high CO2 and warming [11]. On the other hand, it has been clearly recognized that warming is the main driver for the evolution of the marine phytoplankton, rather than high CO2, and led to the marked divergence of thermal reaction norms, and the evolution driven by warming is also found to be constrained by high CO2 [12, 13]. Furthermore, even in multidriver environments, most of the evolutionary changes in population growth rates were explained by some dominant drivers such as high CO2 and elevated temperature [14]. Nevertheless, the molecular basis for this interaction on an evolutionary scale is virtually unknown.

Here, we address this fundamental knowledge gap by carrying out a 2-year selection experiment with the model marine diatom Phaeodactylum tricornutum. We grew three independent replicate cell lines of P. tricornutum for 2 years under each of the following conditions: control (CO2: 400 µatm; temperature: 15 °C), high CO2 (CO2: 1000 µatm; temperature: 15 °C), ocean warming (CO2: 400 µatm; temperature: 20 °C), and high CO2 combined with warming (CO2: 1000 µatm; temperature: 20 °C). After the 2 years of selection, we re-sequenced the genomes and performed a transcriptomics analysis of the populations under those four selection regimes. By applying multi-omics approaches together with various physiological measurements, we aimed to investigate the molecular aspects and/or metabolic pathways that facilitate or constrain the adaptation to high CO2 and/or warming. Therefore, our study may help us in understanding the adaptation capacity of the phytoplankton and in predicting their fate in an acidifying and warming ocean.

Results

Genome map variation

To investigate the underlying mechanisms for the adaptation of P. tricornutum to long-term high CO2 and/or warming selection, we re-sequenced the genomes of the populations in ambient, high CO2, warming, and high CO2+ warming conditions after ~2 years of selection. We identified candidate single-nucleotide polymorphism (SNPs) in populations evolved under high CO2, warming, and high CO2+ warming conditions, compared to those selected in control conditions, that might have contributed to the adaptation to high CO2 and/or warming selection. We identified ~1.9 × 105 SNPs genome wide in P. tricornutum (Table 1). We found that P. tricornutum populations under four selection regimes share most of their SNPs (~1.9 × 105 SNPs, 98.6–99.7% of total) (Supplementary Fig. S1). We also detected pronounced changes in allele frequency of SNPs between the control group with any of the three evolved groups (Supplementary Fig. S2). A total of 1282, 1387, and 1327 SNPs showed distinct allelic distribution in P. tricornutum populations under high CO2, warming, and high CO2+ warming selection, respectively, when compared with that in the control condition (p < 0.05, Fisher test) (Supplementary Fig. S2). Of the ~1.9 × 105 SNPs, most were located in exonic regions (46%), with 9.7% located in coding sequences (Table 1). Among coding regions, there were ~4.6 × 104 synonymous and ~4 × 104 nonsynonymous SNPs with subtle differences among all the four groups (Table 1). An average of ~1.5 × 104 insertions and deletions (indels) were identified in all four groups. Only 18% of indels were located in coding regions, of which 38 and 59% resulted in frameshift and non-frameshift mutations, respectively (Table 1).

Table 1 Summary of the numbers of genomic variants in Phaeodactylum tricornutum populations that were selected by long-term high CO2 and/or warming conditions.
Full size table

Genetic diversity loss after long-term high CO2 and/or warming adaptation

To detect the changes in genetic diversity of P. tricornutum populations after long-term high CO2 and/or warming adaptation, we first calculated the nucleotide diversity π in control and in the other three new evolved conditions. The π of P. tricornutum populations in control conditions (3.230 × 10−3) was higher than that in the three new evolved conditions, with the lowest π in the warming-alone condition (2.964 × 10−3), indicating that a large amount of genetic diversity in P. tricornutum populations had been lost during long-term high CO2 and/or warming selection (Supplementary Fig. S3). Furthermore, a lower π in warming-alone-selected populations than that in high CO2-selected populations suggested that warming is a stronger driving force for the evolution of P. tricornutum, than high CO2. However, this strong selection by warming alone was counteracted by high CO2, showing a weaker selection when high CO2 and warming acted together (3.05 × 10−3), suggesting an antagonistic interaction between high CO2 and warming (Supplementary Fig. S3). We then calculated the population divergence (FST, a genetic measurement of population differential, the ratio of FST in new evolved populations to that under control conditions) to explore the population divergence among the four treatment groups. This analysis indicated a large population divergence between the control group with any of the three evolved groups (high CO2: 0.0265; warming: 0.0304; high CO2+ warming: 0.0253) (Supplementary Fig. S3).

Selection regions during the long-term high CO2 and/or warming adaptation

These two genetic summary statistics, the π ratio and the FST of P. tricornutum populations in control conditions to that in high CO2, warming alone and two together, as described above, were applied to detect the regions with reduced genetic diversity (π) and/or elevated genetic differentiation (FST) in P. tricornutum populations after long-term high CO2 and/or warming adaptation. Specially, we identified a total length of 1.98, 1.63, and 2.70 Mb regions in the P. tricornutum genome with extremely low levels of genetic diversity in the top 5% π ratio values after long-term selection under high CO2 (π ratio ≥1.381), warming (π ratio ≥1.537), and the two together (π ratio ≥1.295) (Fig. 1a). In parallel with these data, a total length of 1.93, 1.53, and 1.68 Mb regions in the P. tricornutum genome with increased levels of genetic differentiation in the top 5% FST values were found associated with selection under long-term high CO2 (FST ≥ 0.063), warming alone (FST ≥ 0.085), and the two factors together (FST ≥ 0.063), respectively (Fig. 1a). Only autosomal regions in the top 5% for both π ratio and FST values were considered candidate regions (i.e., candidate selection regions, CSRs) for adaptation to the new evolved conditions. By applying these thresholds, 434, 450, and 516 genes were identified in autosomal CSRs corresponding to the long-term high CO2, warming, and the combination of high CO2 with warming selection, respectively (Fig. 1d, e and Supplementary Data S1–S3).

Fig. 1: Genome-wide screening of selective sweeps.
figure 1

a–c Genome-wide distribution of selective sweeps identified by FST and π ratio (sliding window = 100 kb, step = 10 kb) during the evolution of Phaeodactylum tricornutum to high CO2 (a), warming alone (b), and its combination with high CO2 (c). Dashed lines with each panel (a–c) indicate the threshold of the whole-genome level (top 5% of FST or π ratio). Only autosomal regions in the top 5% for both π ratio and FST values were considered candidate regions for adaptation to the new evolved conditions. d–f Number of genes encompassed in the selective sweeps during the evolution of Phaeodactylum tricornutum to high CO2 (d), warming alone (e), and its combination with high CO2 (f). The x-axis represents a total of 33 chromosomes in Phaeodactylum tricornutum genome.

Full size image

Evolution of metabolic pathways driven by high CO2 and/or warming adaptation

We hypothesize that the adaptations to high CO2 and/or warming could have been achieved through the simultaneous expressions of functionally related genes in CSRs. To test this hypothesis, we performed enrichment analysis on gene ontology (GO) terms. With this analysis, we found that the GO terms and the majority of the associated genes (top 3) within the adaptations of P. tricornutum populations to high CO2, warming alone, and warming combined with high CO2 were quite similar. In general, we identified 35–36 GO terms, and most of the genes (top 3) were related to “metabolic process (high CO2: 110 genes; warming: 112 genes; high CO2+ warming: 130 genes)”, “catalytic activity (high CO2: 97 genes; warming: 99 genes; high CO2+ warming: 123 genes)”, “cellular process/binding (high CO2: 89 genes; warming: 96 genes; high CO2+ warming: 112 genes)” (Supplementary Figs. S3–S5 and Supplementary Data S1–S3). In addition, genes involved in cellular component organization or biogenesis, signaling, cell junction, and protein binding were also identified (Supplementary Fig. S4 and Supplementary Data S1–S3).

To identify how P. tricornutum populations adapt to high CO2, warming, and warming combined with high CO2 conditions via altering metabolic pathways, the genes in CSRs were mapped to the known metabolic pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) (Supplementary Figs. S5–S7). For the long-term high CO2 selection, our results showed that the gene for the triose-phosphate transporter (TPT), which catalyzes the exchange of 3-phosphoglycerate (3-PGA), triose-phosphate, and inorganic phosphate across the plastid envelope, and plays vital roles in photosynthesis, was affected (Fig. 2a and Supplementary Table S4). In parallel with TPT, the gene for phosphoglycolate phosphatase (PGLYCP), which catalyzes the conversion from 2-phosphoglycolate to glycolate in photorespiration, a pathway that competes with photosynthetic carbon fixation, was also affected (Fig. 2a). Glycolate is oxidized in the peroxisome and is one of the intermediates in the glyoxylate pathway (Fig. 2a). Furthermore, the gene for isocitrate lyase (ICL), which catalyzes the reaction from citrate to isocitrate in the glyoxylate pathway in the peroxisome was also affected by the long-term high CO2 adaptation (Fig. 2a). In addition, we identified the gene for methylmalonate-semialdehyde dehydrogenase (MMSDH), which functions in branched-chain amino acid (BCAA) degradation, as being affected after long-term high CO2 adaptation (Fig. 2a). As the glyoxylate pathway shares some of the intermediates (e.g., citrate, isocitrate) with the TCA cycle and succinyl-CoA, and the product of BCAA degradation is also one of the important intermediates in the TCA cycle, it is expected that TCA would be altered by long-term high CO2 selection. In line with our expectations, we found that the gene for 2-oxoglutarate dehydrogenase (OGDC), which functions in the catalysis of α-ketoglutarate to succinyl-CoA in TCA, was affected (Fig. 2a). Along with TCA cycle enzymes, the genes for phosphoglycerate kinase (PGK) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH), which function in glycolysis, another key metabolic pathway for energy-production along with the TCA, were affected by high CO2 (Fig. 2a). Glycolysis breaks down glucose and produces (intermediate) metabolites (e.g., acetyl-CoA) to participate in the TCA cycle and fatty acid biosynthesis (Fig. 2a). In line with this, our results showed that the gene for 3-oxoacyl-acyl-carrier-protein reductase (FABG), which is responsible for the extension of fatty acids from malonyl-ACP in fatty acid metabolism, and the gene for triacylglycerol lipase (ATGL), which is responsible for the initial reaction of triacylglycerol (TAG) breakdown, were affected after long-term high CO2 adaptation (Fig. 2a). The observation that genes for enzymes that function in key energy-production pathways (e.g., TCA cycle, glycolysis) were affected by high CO2 adaptation suggests that energy metabolism evolved rapidly to the new high CO2 conditions. This finding was further supported by the high CO2 selections on the genes for (mitochondrial) ATP synthase (ATPA), adenylate kinase (ADK), and fucoxanthin chlorophyll a/c protein (FCPs) (Fig. 2a). The enzymes encoded by the former two genes (i.e., ATPA and ADK) catalyze the reversible conversion between ATP and ADP, and the latter protein has been shown to play a critical role in light harvesting in diatoms.

Fig. 2: Evolution of metabolic pathways.
figure 2

a Metabolic pathways including photosynthesis, photorespiration, oxidative pentose phosphate pathway, triacylglycerol (TAG) accumulation, glycolysis, fatty acid metabolism, glyoxylate pathway, ATP metabolism, oxidative phosphorylation, and diel variability, driven by genes selected by high CO2 (blue), warming (blue + cyan) (i.e., gene names with blue or cyan were affected) and its combination with high CO2 (blue + cyan + purple) (i.e., gene names with blue or cyan or purple were affected). Cyan symbols with asterisk indicate that the gene is exclusively affected by warming. b–d Heatmap displays the log2-fold change in gene expressions between control and high CO2 (b) or warming (c) or its combination with high CO2 d selected populations as presented in a. It should be noted that the second half of the glycolysis reaction chain (from G-3-P to pyruvate) has been duplicated, and the respective isoenzymes can also be found in mitochondria and plastid. However, here we only show the glycolysis in cytosol to avoid any redundancy. PGLYCP phosphoglycolate phosphatase, TPT triose-phosphate transporter, FCPA/FCPE fucoxanthin chlorophyll a/c protein, XI xylose isomerase, TKT transketolase, ATGL triacylglycerol lipase, HAD 3-hydroxyacyl-[acyl-carrier protein] dehydrase, FABB 3-oxoacyl-[acyl-carrier-protein] synthase, FABG 3-oxoacyl-acyl-carrier-protein reductase, OGDC 2-oxoglutarate dehydrogenase, ATPA ATP synthase, ADK adenylate kinase, PCC propionyl-CoA carboxylase, MMSDH methylmalonate-semialdehyde dehydrogenase, PGK phosphoglycerate kinase, GAPDH glyceraldehyde-3-phosphate dehydrogenase, CK casein kinase, G-3-P glyceraldehyde-3-phosphate.

Full size image

In terms of evolution to warming, we found that the genes that were mapped to the known metabolic pathways in the KEGG, affected by long-term high CO2 as described above, were also represented in the long-term warming-selected populations (Fig. 2a and Supplementary Data S4). In addition to this, we identified the gene for transketolase (TKT) involved in the oxidative pentose phosphate pathway (OPPP), which produces NADPH and pentose phosphate for the biosynthesis of amino acids, nucleotides, and fatty acids by decarboxylating glucose-6-phosphate, was also affected by warming (Fig. 2a). In line with this finding, we found that gene for NADH dehydrogenase (NDH), which accepts electrons from NADH+ and H+ derived from other metabolic pathways (e.g., the OPPP) to create an electrochemical gradient across the inner mitochondrial membrane as part of the respiratory electron transport in oxidative phosphorylation, was affected by long-term warming (Fig. 2a).

We then investigated how the genes might evolve to the long-term high CO2 combined with warming selection. It was found that the gene for NADP dehydrogenase (NDH), affected under the warming-alone treatment, was not affected under warming combined with high CO2 conditions, suggesting an antagonistic interaction between high CO2 and warming selection (Fig. 2a). In addition to the warming-alone selected genes that are involved in the metabolic pathways, as shown above, we identified two genes: those for 3-hydroxyacyl-[acyl-carrier protein] dehydrase (HAD) and 3-oxoacyl-[acyl-carrier-protein] synthase (FABB), which function in fatty acid degradation and biosynthesis, as being affected by the warming combined with high CO2 condition (Fig. 2a). The gene for propionyl-CoA carboxylase (PCC), which functions in BCAA degradation, and the gene for xylose isomerase (XI), which is involved in xylulose-5-P (an intermediate of OPPP) metabolism were also affected (Fig. 2a).

Genomics and transcripts response interaction

To test whether the genes involved in the metabolic pathways at the genomic level, as described above, were expressed differentially at the mRNA level, we then further performed a transcriptomics analysis. Although the results showed that the relationships between the changes in allele frequency of SNPs in exons and changes in gene expressions were weak at the whole genomic level in all the three evolved populations (Supplementary Fig. S8), we found that there was an overlap (10–24%) between the genes selected by high CO2 and/or warming conditions and the genes differentially expressed after high CO2 and/or warming conditions (Supplementary Data S5–S7). We found that 6 and 36 genes out of 430 affected genes in high CO2-alone populations were significantly up- and downregulated, respectively (Supplementary Data S5). For warming-selected populations, 14 and 66 genes out of 450 affected genes were significantly up- and downregulated, respectively (Supplementary Data S6), while for gene expression in the populations selected under warming combined with high CO2, we found 36 genes were upregulated and 86 genes were downregulated out of the total of 516 affected genes (Supplementary Data S7). For the genes involved in the metabolic pathways or biological processes described above, our results showed that while selection under ocean warming alone significantly decreased the gene expression of ICL and NDH, there was no difference in the populations selected by warming combined with high CO2 (Fig. 2b–d and Supplementary Data S4). In addition, the insignificantly expressed genes for MMSDH, ATGL, and PGLYCP found in warming-alone selected populations were significantly downregulated when warming selection was combined with high CO2 (Fig. 2b–d and Supplementary Data S4). Furthermore, warming-alone selection led to larger downregulation of gene expression of PGK, GAPDH, and ATPA compared to that of high CO2-alone selection, and these strong inhibitory effects were counteracted by high CO2 selection (i.e., an antagonistic interaction), showing relative smaller downregulation in the populations selected under warming combined with high CO2 (Fig. 2b–d and Supplementary Data S4). While the gene expression of TPT was upregulated in warming-alone selected populations, it was downregulated in high CO2-alone or in its combination with warming selected populations (Fig. 2b–d and Supplementary Data S4).

Evolution of metabolic traits driven by high CO2 and/or warming adaptation

We hypothesize that the potential changes to metabolic pathways as a result of genomic variants and gene expressions at the transcriptomics level are also linked to some key metabolic traits at the physiological level. Therefore, we determined the metabolic traits of photosynthesis and respiration, two of the most important metabolic traits that may affect biogeochemical cycles, for the populations in the four selection regimes. Consistent with the results at the genomics level, a stronger selection by warming resulted in significantly higher rates of gross photosynthesis (24%) and respiration (20%) (linear mixed-effects model [LME], p < 0.05) (Fig. 3a, b and Supplementary Tables S1 and S2). As nearly the same enhancement of rates of gross photosynthesis and respiration were found under warming selection, no significant differences in carbon usage efficiency (CUE) were detected (LME, both p > 0.05) (Fig. 3c and Supplementary Tables S1 and S2). These large enhancements in metabolic traits under warming selection were expected to be linked to genomics and transcriptomics. However, expression of the gene for NDH, which was specifically affected by warming and functions in transferring electrons to ubiquinone in the respiratory electron transport chain, was significantly downregulated (log2(fold change) = 1.35) (Fig. 2c and Supplementary Data S4).

Fig. 3: Evolution of metabolic traits.
figure 3

a Gross photosynthetic rates (pmol O2 cell−1 h−1), b respiration rates (pmol O2 cell−1 h−1), and c the carbon-use efficiency (CUE) in Phaeodactylum tricornutum populations after 2 years high CO2 and/or warming selection. The different letters indicate significant differences between treatments tested by the linear mixed-effects model (LME, details see in Materials and methods).

Full size image

High CO2 alone or combined with warming adaptation alters the diel variability

Besides the genes involved in the processes described above, the GO and KEGG analyses also pinpointed that one gene, namely that for casein kinase (CK) involved in circadian rhythms [15], was affected under all the three evolved conditions (Fig. 2a and Supplementary Data S1–S3). Consistent with the findings at the genomics level, the transcriptomics data showed that there were in total five differentially expressed circadian rhythm genes in all the three new evolved conditions when compared with expression in control conditions (Fig. 4a–c and Supplementary Data S5–S7). Specially, we found that three of these genes (CRY2, CSK2B, and CK2B-1) were significantly upregulated (log2(fold change): ~1.40–1.55, all p < 0.05) in long-term high CO2 selected populations, and only one gene, namely WDR61, was significantly downregulated (log2(fold change) = −0.9, p = 0.047) in the warming-selected populations (Fig. 4a, b and Supplementary Table S10). However, when warming and high CO2 acted together, the three upregulated genes (CRY2, CSK2B, and CK2B-1, all p < 0.05) in high CO2-alone selected populations were relatively downregulated (log2(fold change): ~1.0–1.3, all p < 0.05) (i.e., an antagonistic interaction), and the gene WDR61 was further downregulated (log2(fold change) = −1.29) (Fig. 4c and Supplementary Data S8). Given the functions of these four genes, these two lines of evidence from genomics and transcriptomics imply that the long-term high CO2 and/or warming adaptation might alter the diel variability of P. tricornutum at the physiological level. To test this hypothesis, we measured the photochemical traits of the rapid light curve (RLC) and maximal photochemical quantum yield of PSII (Fv/Fm) for the populations under the 4 selection regimes at 2-h intervals over a 24-h diel cycle. The results showed that long-term high CO2 selection altered the diel cycle of Fv/Fm, and that of relative maximal relative electron transport rate (rETRmax), saturating light intensity Ik, and apparent photosynthetic efficiency α were significantly decreased by warming selection (Fig. 4d. GAM model and Supplementary Table S3). The model also indicated a significant interaction of high CO2 and warming selection on the apparent photosynthetic efficiency α (Fig. 4d. GAM model and Supplementary Table S3). It is worth noting that a shift in saturating light level in populations under warming and its combination with high CO2 selection may exert influences on other photochemical responses, and lack of reciprocal transplant assay may not allow us to test whether the photochemical changes in diel variability were a direct response to warming or its combination with high CO2 selection. However, even if the photochemical changes in diel variability observed are directly due to the shifts in light requirement, they are originally driven by warming or its combination with high CO2 selection as the populations under the four selection regimes were grown under the same conditions (e.g., light, nutrients) except for CO2 and temperature. Therefore, these physiological results together with the multi-omics data support the contention that high CO2 and its combination with warming selection altered the diel variability of the P. tricornutum populations.

Fig. 4: The alterations of diel variability in Phaeodactylum tricornutum.
figure 4

a Heatmap displays the log2-fold change in gene expressions between control and high CO2 (a) or warming (b) or its combination with high CO2 c selected populations involved in diel variability. d Diel changes in photochemical traits of maximal photochemical quantum yield of PSII (Fv/Fm), relative maximal relative electron transport rate (rETRmax), saturating light intensity (Ik), and the apparent photosynthetic efficiency (α) in the cells under different selection conditions. Data are the mean ± SE (n = 3).

Full size image

Discussion

While a large body of literature has reported that high CO2 and warming could act additively, antagonistically, or synergistically on various marine organisms at short-term acclimatory and long-term adaptive scales [12, 13, 16, 17], the molecular basis for such interactions in an evolutionary context has rarely been reported. In the present study, the lowest genetic diversity and the highest genetic differentiation in the populations after long-term warming selection showed that ocean warming is a stronger driving force for the evolution of P. tricornutum than high CO2. However, a recovery of genetic diversity and genetic differentiation in the populations exposed to long-term high CO2 and warming selection together, suggested that the evolution driven by warming can be constrained by high CO2. We then further elucidated the underlying molecular basis for this constraint and found that high CO2 selection altered the expression of warming-affected genes involved in some key metabolic pathways or biological processes, such as the glyoxylate pathway, amino acid, fatty acid metabolism, and diel variability. Diatoms are reported to sense a change in external CO2 levels via cAMP signaling and coordinated gene expression [18,19,20], in which the enhanced production of cAMP induces the CO2 concentrating mechanism (CCM)/photorespiration cluster transcription factors to bind to the TGACGT motif and then inhibit the transcription of CCM/photorespiration cluster genes [20]. Given the alterations of high CO2 selection on the warming-affected genes involved in the metabolic pathways or biological processes described above, it is an interesting open question whether the cAMP signaling system is responsible for these alterations. While it is not surprising that warming, rather than high CO2, is the main driver for the observed evolution as most enzymes are temperature dependent [11, 21], and that the evolution driven by warming can be constrained by high CO2 as demonstrated in our previous study [13], here we provide a holistic understanding of the interactive effects of high CO2 and warming on phytoplankton in an evolutionary context by integrating genomics, transcriptomics. and physiological data. Marine organisms are exposed to multiple environmental changes, including high CO2, warming, and deoxygenation [3, 22], thus identifying the evolutionary responses of marine organisms to multifaced marine environmental changes and distinguishing the underlying molecular mechanisms for the interactions among multiple drivers are challenging [17]. Hence, our results shed new light on the evolutionary responses of marine phytoplankton to multiple environmental changes and provide new insight into the molecular basis underpinning interactions among those multiple drivers.

Our results demonstrate the potential association of the ability of marine diatoms to adapt to high CO2 and/or warming conditions with genomic variants and differential gene expression, which can lead to substantial changes in some key metabolic traits. A larger genetic diversity loss and a higher genetic differentiation in the populations adapted to warming than that to high CO2 led to a larger enhancement in metabolic traits of photosynthesis and respiration. These findings suggest that the observed changes in genomic variants and gene expressions were also linked to the changes in metabolic traits, as demonstrated in some of the previous studies [8, 23, 24]. However, this consistency between studies was not always true for growth, a common fitness proxy for experimental evolution studies [7]. For instance, while hundreds of genes were affected by high CO2 and/or warming conditions, there were no significant changes in the functional trait of growth at the end of the long-term selection in the present study (Supplementary Fig. S9). Because phenotype, not the genotype, is the direct target of selection, with a mutant’s fitness as the ultimate target [25], evolutionary trait changes were generally measured and defined using trait changes that do not revert when the population is placed back in the control environment (i.e., a reciprocal transplant) [26]. Hence, the lack of reciprocal transplant assays in the present study may not allow us to quantify the evolutionary trait changes. However, one key consideration is that trait changes do not always correlate with the gene changes at the genomic level as found with the inconsistency of growth as observed in the present study and other studies [27]. For instance, a larger downregulation of gene expression of PGK, GAPDH, and ATPA in warming-alone selected populations than that in high CO2-alone selected populations are expected to cause a lower respiration rate [28], which in fact is higher in the present study, coupling with a higher photosynthesis rate. A higher photosynthesis rate in warming-alone-selected, and its combination with high CO2-selected, populations is expected to be linked to enhanced gene expression of TPT, which functions in the exchange of 3-PGA, triose-phosphate, and inorganic phosphate across the plastid envelope, and plays vital roles in photosynthesis, in populations under these two corresponding selections regimes. However, this gene was only upregulated in the warming-alone selected populations. Therefore, our study suggests that the adaptations of marine phytoplankton to ongoing marine environmental changes can be achieved through simultaneous changes in genomic variants, coordinated gene expression, and corresponding metabolic traits, but that care needs to be taken when interpreting the intrinsic correlations among these changes.

Most ecosystem models of the response of marine phytoplankton to global change do not consider their adaptive capacities, thus predictions of phytoplankton distribution and likely changes in global ocean primary production are governed by our understanding of short-term responses and acclimation to the environmental conditions associated with global change [11, 29]. However, recent studies have presented evidence for the evolution of marine phytoplankton in adaptive responses to marine environmental change [8, 9, 28]. Therefore, the outcomes of long-term evolutionary studies can be used to predict how traits that affect biogeochemical cycles (such as photosynthesis and respiration) and marine food webs (such as cell size) may evolve in future oceans [7, 30, 31]. The long-term evolutionary response traits of such as photosynthesis and respiration in relation to high CO2 and/or warming in the present study can therefore be integrated into models of ocean biogeochemistry and improve projections of marine primary production for the coming decades. Furthermore, the trait-based genomics and transcriptomics data presented in our study might also be used to track ocean biogeochemistry as proposed by Coles et al. [32]. Taken together, our study suggests that the integrative approach by examining the evolutionary responses of marine phytoplankton at multi-omics and physiological levels may provide more accurate projections of ocean biogeochemistry over the twenty-first century.

Although novel mutations have been shown to be the genetic source of adaptation as demonstrated in various empirical studies [33, 34], the importance of standing genetic variation in the adaptation of species to changing environments has also been highlighted in some previous work [35, 36]. However, as our experiments were conducted with a single clone of P. tricornutum, we were unlikely to identify the standing genetic variation which can contribute substantially to the adaptation of marine diatoms with highly diverse meta-populations [37]. Furthermore, it has been reported that P. tricornutum has strong genetic subdivisions of the populations sampled at broad geospatial and temporal scales, and showed compelling signatures of genetic and functional convergence [38]. Consequently, our study with a single clone of P. tricornutum may be a conservative estimate of the evolutionary potential of this species in coping with the ongoing marine environmental changes. In addition, while no samples were taken at the beginning of long-term selection experiments for genome sequencing, we would be not able to identify the mutations arising de novo after the selection initialized, and becoming fixed at the end of the selection experiment for the populations under each selection regime. Although this type of genetic adaptation has been suggested to contribute to a quite low proportion of the total variations [8], it is still worth measuring to get a complete picture of the genetic changes that drive adaptation to a changing environment. Different phytoplankton species between- and within-taxonomic groups exhibited various adaptive capacities in response to ongoing marine environmental changes [39, 40]. Therefore, we propose that experiments on meta-populations across diverse phytoplankton taxonomic groups are highly warranted to improve estimations of the evolutionary potentials of marine phytoplankton to global change.

Materials and methods

Culture conditions

A monoclonal culture of P. tricornutum Bohlin bac-2, obtained from the Seaweed Culture Collection Center, Institute of Oceanology, Chinese Academy of Sciences (No: MASCC-0025, publicly available upon request), was maintained in half-strength Guillard’s “F” solution [41], containing 100 μmol L−1 silicate, 100 μmol L−1 nitrate, and 6.25 μmol L−1 phosphate under a photon flux of 100 μmol photons m−2 s−1 with a light: dark cycle of 12 h:12 h (HP1000G-D, Ruihua) in 15 °C plant growth chambers (HP1000G-D, Ruihua). The cells were synchronized by prolonged darkness (48 h) before the long-term selection experiments [42]. For the long-term selection experiments, the single-clone cultures were diluted to separate biological replicates and were grown in 500 mL Nalgene bottles at control (CO2: 400 µatm; temperature: 15 °C), high CO2 (CO2: 1000 µatm; temperature: 15 °C), ocean warming (CO2: 400 µatm; temperature: 20 °C), and high CO2 combined with warming conditions (CO2: 1000 µatm; temperature: 20 °C). The choice of 20 °C and 1000 µatm as elevated temperature and pCO2, represent the projected values at the end of this century according to the IPCC high emission scenario RCP SSP5-8.5 [43]. Although some previous studies have reported that 20 °C is close to the optimal growth temperature of this species [44, 45], our previous work demonstrated that 20 °C is a stressor for the strain used in the present study due to strain-specific thermal reaction norms of diatoms [46, 47]. The CO2 level of 400 µatm was reached by pre-aerating the medium with the ambient outdoor air, while the CO2 level of 1000 µatm was achieved within a plant growth chamber (HP1000G-D, Ruihua), in which the target CO2 levels were obtained by mixing air and pure CO2 gas, and the CO2 partial pressure was continuously monitored and maintained at 1000 ± 50 µatm. Three independent replicate semi-continuous batch cultures were run for ~2 years under each of the four selection regimes (i.e., control, high CO2, warming, high CO2+ warming). The initial cell concentration was 50 cells mL−1 and the cultures were partially renewed with fresh medium equilibrated with the corresponding target CO2 levels every 5–7 days to restore the cell density to the initial level. The cell densities were thereby maintained within a range of ~4.0 × 104–6.0 × 105 cells mL−1 at the time of dilution. Nutrients were not limiting as the cell abundances achieved at the end of the batch cycles were far from those expected at the stationary phase (~2.0 × 106 cells mL−1). To avoid any gas exchange, all the cultures were kept in tightly closed polycarbonate bottles that were filled with culture medium. Although we were unable to determine the carbonate chemistry parameters weekly for a 2-year experiment, our pilot experiments indicated that this semi-continuous culture approach is reliable [28]. By multiplying the number of generations for the cells under control or high CO2 selection after ~400 days of selection as reported in our previous study [28] and 1.825 (the conversion factor from 400 days to 2 years), we estimated the approximate number of generations for this diatom species after 2 years selection in control or high CO2 condition. Here, we show that after the 2 years selection period, P. tricornutum had grown for approximately 1615 and 1611 generations under control and high CO2 treatments, respectively. And P. tricornutum populations had grown for approximately 1639 and 1644 generations after 2 years of selection under warming alone or warming combined with high CO2. The populations under 4 selection regimes showed fluctuations in specific growth rate (e.g., an increase in growth rate from ~day 273 to 300 and a decrease in growth rate from ~day 300 to 320) (Supplementary Fig. S9).

DNA extraction and whole-genome re-sequencing

At the end of the 2-year long-term selection experiments, the triplicate (~109 cells for each replicate) samples from control, high CO2, warming, or high CO2+ warming selection regimes at the same time point (i.e., same cell cycle stage) in the middle of light phase were pooled by selection regime and then used for DNA extraction. Approximately 3 µg genomic DNA was used to construct paired-end libraries with an insert size of 500 bp using a Paired-End DNA Sample Prep kit (Illumina Inc., San Diego, CA, USA). These libraries were sequenced using a NovaSeq 6000 (Illumina Inc., San Diego, CA, USA) NGS platform at Gene Denovo Biotechnology company (Guangzhou, China). Quality trimming is an essential step to generate high confidence in variant calling. Raw reads were processed to get high-quality clean reads according to the following stringent filtering standards of fastp (version 0.18.0) [48]: (1) removing reads containing adapters; (2) removing reads containing more than 10% of unknown nucleotides (N); (3) removing low quality reads containing more than 50% of low quality (Q-value ≤20). The resulting reads were then mapped against the existing P. tricornutum genome [49]. The average sequencing depth for populations under four selection regimes is 60× to 100×, and this is sufficient to cover the whole genome of this species (Supplementary Fig. S10).

Variant identification and annotation

We identified SNPs and InDels in Burrows–Wheeler Aligner of genomic sequence reads against the reference genome sequence with the settings “mem 4-k 32 -M”. Here, -k is the minimum seed length, and -M is an option used to mark shorter split alignment hits as secondary alignments [50]. Variant calling was performed for all samples using the GATK’s Unified Genotyper. SNPs and InDels were filtered using GATK’s Variant Filtration with proper standards (-Window 4, -filter “QD < 2.0 || FS > 60.0 || MQ < 40.0 “, -G_filter “GQ < 20”) and those exhibiting segregation distortion or sequencing errors were discarded. The physical position of each SNP was determined by the software tool ANNOVAR [51].

Selective sweep analysis

A sliding-window approach (100-kb windows sliding in 10-kb steps) was applied to quantify polymorphism levels (π, nucleotide diversity as a measure of variability) and genetic differentiation (FST) between P. tricornutum populations in control and the three new evolved conditions. To detect regions with significant signatures of a selective sweep, we considered the distribution of the π ratios (πcontrol/πevolved) and FST values. Only autosomal regions in the top 5% for both π ratio and FST values were regarded as candidate regions for adaptation to the new evolved conditions, and genes in these regions were considered candidate genes.

Gene Ontology enrichment and KEGG analysis

All the candidate genes were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/). The significantly enriched GO terms were defined by a hypergeometric test. The GO project classifies genes into three categories based on their annotation: molecular function, biological process, and cellular component. Following the annotation, we further performed a pathway-based analysis to identify the significantly enriched metabolic pathways or signal transduction pathways for the candidate genes in the KEGG (http://www.genome.jp/kegg/).

Transcriptome analysis

In parallel with the genomic analysis, three biological replicate samples were also collected for transcriptome analysis at the end of the 2-year long-term selection experiments. Briefly, 1L P. tricornutum cultures at the relevant selection regime were collected at the same time point (i.e., same cell cycle stage) in the middle of the photoperiod by centrifugation (8000 g, 10 min), and the samples were then flash-frozen in liquid nitrogen and stored at −80 °C until further analysis. Total RNA was extracted using a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase-free agarose gel electrophoresis. After the total RNA extraction, the eukaryotic mRNA was enriched by Oligo (dT) beads, while the prokaryotic mRNA was enriched by removing ribosomal RNA (rRNA) by a Ribo-Zero Magnetic Kit (Epicentre, Madison, WI, USA). The enriched mRNA was then fragmented into short fragments using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP, and buffer. The cDNA fragments were then purified with a QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end-repaired, poly (A) added, and ligated to Illumina sequencing adapters. The ligation products were size-selected by agarose gel electrophoresis, PCR amplified, and sequenced using NovaSeq 6000 by Gene Denovo Biotechnology company (Guangzhou, China).

To get high-quality clean reads, the reads obtained from the sequencing machines were filtered by fastp (version 0.18.0) [48]. The short reads alignment tool Bowtie2 (version 2.2.8) was used for mapping reads to a rRNA database [52]. The rRNA mapped reads then were removed, and the remaining clean reads were then mapped against the P. tricornutum genome [49] using HISAT (version 2.2.1) [53]. The mapped reads of each sample were assembled using StringTie (version 1.3.1) in a reference-based approach [54, 55]. And then FPKM (fragment per kilobase of transcript per million mapped reads) values were calculated for each gene to quantify its expression abundance and variation. RNA differential expression analysis was performed by DESeq2 software between two different groups [56]. The genes with values of the parameters false discovery rate below 0.05 and absolute fold change ≥2 were considered differentially expressed genes.

Growth rate measurement

The growth of P. tricornutum populations in the four regimes was determined at the end of the 2-year selection experiments. For the determination of growth, cell number was quantified by examining the samples under an optical microscope (Optec B203LEDR, Optec, China) by a hemocytometer. The specific growth rate, μ (d−1), was calculated as follows:

$$\mu = \frac{{{{{{{{{\mathrm{ln}}}}}}}}({{C}}_1/{{C}}_0)}}{{({{t}}_1 - {{t}}_0)}}$$

where C1 and C0 are the cell densities at times t1 and t0, respectively. The number of generations per transfer (g) is equivalent to the number of doublings and was calculated as follows:

$${{g}} = \frac{{{{t}}_1 - {{t}}_0}}{{{{{{{{{\mathrm{ln}}}}}}}}\left( 2 \right)/{{{{{{{\mathrm{\mu }}}}}}}}}}$$

where t1 − t0 is the time interval of the transfer (d), ln (2)/μ is the doubling time (d), and μ is the specific growth rate (d−1). The samples taken for specific growth rates were all collected in the middle of the light phase.

Determination of photosynthesis and respiration

At the end of the long-term selection experiments, populations under the four selection regimes were harvested in the middle of the light period and concentrated by filtration. The concentrated cells then were resuspended into 20 mmol L−1 Tris-buffered medium whose pH was adjusted, using freshly prepared hydrochloric acid and sodium hydroxide, to their corresponding culture medium values. Photosynthetic (and respiration) rates were thus measured at CO2 and bicarbonate levels identical to those for growth. Net photosynthesis (P) was measured through oxygen evolution in the light (100 µmol photons m−2 s−1, equal to the growth light level) using a Clark-type oxygen electrode (YSI 5300A; Yellow Springs Instrument Co., Inc., OH., USA). Respiration (R) was determined by measuring O2 uptake in darkness for 10 min after the samples were dark-adapted for ~15 min. The rates of P and R at the temperature used for growth were normalized by dividing by the cell numbers measured in each aliquot (pmol O2 cell−1 h−1).

Gross photosynthesis was estimated as Pgross = R + P. We also calculated the CUE, which represents the fraction of fixed carbon that is available for allocation to growth and can be estimated from the rates of photosynthesis (P) and respiration (R) [57] as follows:

$${{{{{{{\mathrm{CUE}}}}}}}} = 1 - {{R}}/{{P}}_{{{{{{{{\mathrm{gross}}}}}}}}}$$

CUE was estimated for each independent replicate from rates of R and P measured at their selection temperature and DIC.

Quantifying diel rhythm of photochemical responses

To assess the diel rhythm of photochemical responses of P. tricornutum populations in the four selection regimes, we measured the maximum quantum yield of photosystem II (Fv/Fm) and RLCs using a pulse amplitude modulated fluorometer (PAM-2100, Waltz, Effeltrich, Germany) at 2-h intervals over 24 h. For the determination of Fv/Fm, all samples were dark-adapted for 15 min before measurements, and the saturating pulse was applied at 5000 μmol photons m−2 s−1 for 0.8 s. Fv/Fm was calculated as Fv/Fm = (Fm – F0)/Fm, where Fm represents the maximum fluorescence yield after dark adaptation and F0 the minimal chlorophyll fluorescence in the presence of the low intensity modulated measuring light. For the measurements of RLCs, 3 mL samples were incubated at 100 μmol photons m−2 s−1 and 15 or 20 °C according to the temperature in long-term selection regimes for 10 min to avoid induction effects on the photosystems caused by quasi-dark adaptation during manipulation. The RLC is the rETR in response to eight different and increasing actinic photon fluxes (24, 79, 117, 179, 263, 401, 598, and 849 μmol photons m−2 s−1) with a 10 s duration for each increment separated by a 0.8 s saturating pulse (5000 μmol photons m−2 s−1). The rETR (an arbitrary unit) was calculated as follows:

$${{{{{{{\mathrm{rETR}}}}}}}} = {{{{{{{\mathrm{Yield}}}}}}}} \times 0.5 \times {{{{{{{\mathrm{PFD}}}}}}}}$$

where Yield is the effective photochemical quantum yield of photosystem II (PSII) in the light, PFD is the photon flux (μmol photons m−2 s−1), and the factor 0.5 accounts for approximately 50% of all the absorbed energy being allocated to PSII. The RLC was fitted as follows:

$${{{{{{{\mathrm{rETR}}}}}}}} = {{{{{{{\mathrm{rETR}}}}}}}}_{{{{{{{{\mathrm{max}}}}}}}}} \times \tanh \left( {{{{{{{{\mathrm{\alpha }}}}}}}} \times \frac{{{{{{{{{\mathrm{PFD}}}}}}}}}}{{{{{{{{{\mathrm{rETR}}}}}}}}_{{{{{{{{\mathrm{max}}}}}}}}}}}} \right)$$

where rETRmax is the maximal rETR, and α is the apparent photosynthetic efficiency [58]. Ik, the photon flux at which electron transport is initially saturated, was calculated as follows:

$${{I}}_{{k}} = \frac{{{{{{{{{\mathrm{rETR}}}}}}}}_{{{{{{{{\mathrm{max}}}}}}}}}}}{{{{{{{{\mathrm{\alpha }}}}}}}}}$$

Statistical analysis

The interactions between long-term high CO2 selection and warming selection on metabolic traits were analyzed by an LME. In the model, the metabolic trait (e.g., respiration) was considered the dependent variable, high CO2 selection and warming selection were the fixed effects, and biological replicate was treated as a random effect. Model fitting and selection started with the full model and then a series of reduced models with interaction terms and main effects removed. For multi-model selection, we computed small sample-sized corrected Akaike Information Criterion scores (AIC) and then compared between models by calculating delta AICc values and AICc weights using the “dredge” function in “MuMIn” package (see Supplementary Table S5 for all tested models’ AICc). LME was fitted to the data using the “nlme” package and was conducted in R (v. 4.1.1). Multiple comparisons of means were performed on significant effects using the generalized linear hypothesis test (glht) and Tukey’s test in the package “multcomp”.

The photochemical parameters Fv/Fm, rETRmax, α, and Ik were analyzed using a generalized additive mixed-effects model (GAMM) to examine whether the fitness trajectories differed among the four selection regimes. The most complex models included an effect of CO2 level and an effect of temperature and also allowed the shape of the diel rhythm, which was modeled using a cubic regression spline, to vary between selection regimes. The effects of CO2 and temperature on the shape were modeled as fixed effects in the GAMM. Model selection proceeded as described above, which entailed fitting a range of models to the data, starting with the full model and then a series of reduced models with interaction terms and main effects removed. GAMMs were fitted to the data using the “gamm4” package and were conducted in R (v. 4.1.1). We checked for homogeneity of variances and normal distribution and found no major violation of the test assumptions.