Seasonal temperature variation in Zostera marina seedlings under ocean acidification

1 Introduction

As marine angiosperms, seagrasses act as physical barriers to reduce waves impacts, such as turbulence (Granata et al., 2001; Koch et al., 2007). Seagrass beds are recognized as carbon sinks (Fourqurean et al., 2012), effective in accumulating subsurface carbon deposits (Mcleod et al., 2011; Fourqurean et al., 2012) and capturing carbon (Duarte et al., 2005; Barbier et al., 2011), and in addition, seagrass beds can recycle nutrients (Hemminga et al., 1991), provide habitat for many epiphytes and animals (Coll et al., 2011), are nurseries for many fish species, and have high primary productivity (Espino et al., 2011). Seagrass beds can also play a key role in biogeochemical cycling in nearshore waters with these characteristics (Hemminga and Duarte, 2000; Celia and Fernando, 2023). Therefore, seagrass beds with considerable conservation value deserve attention and protection (Barbier et al., 2011).

It has been found that some seagrasses grow rapidly with an annual reproduction pattern (Orth et al., 2006a). The vitality of these seagrasses has a distinct seasonal character (Marbà et al., 1996). Zostera marina, a common seagrass species in the temperate zone of the Northern Hemisphere (Olsen et al., 2016), is a typical dominant species of temperate seagrasses in China, forming a large number of Z. marina meadows along the coast of the Yellow Sea and Bohai Sea in China (Xu et al., 2019). Yan et al. (2023) found that the sea surface temperatures (SST) in the Yellow Sea and Bohai Sea in China over the past 20 years were highly seasonal. It was shown that such seasonal temperature variations affect the growth and productivity of vegetation in marine environments (Sand-Jensen, 1975; Orth and Moore, 1986; Watanabe et al., 2005). And the growth of coastal seagrasses is greatest in summer, with increases in biomass, shoot density and canopy height (Boström et al., 2003). As winter temperatures drop, seagrass leaves virtually disappear, reducing aboveground biomass and eventually leaving only belowground rhizomes to acclimatize to cold temperatures at the lowest of growth (Baden and Pihl, 1984; Burke et al., 1996).

Seasonal changes in the local environment of coastal areas imply not only SST changes but also in other abiotic variables. For example, SST, light, pH, and sediments can regulate seagrass growth (Hemminga and Duarte, 2000; Lee et al., 2005). Correlations exist between environmental variables (Celia and Fernando, 2023). Specifically, there are spatial pattern correlations across latitudinal gradients among seawater temperature, nutrients concentrations, and plants primary productivity (Parmesan and Yohe, 2003; Lemos and Pires, 2004; Tuya et al., 2012; Celia and Fernando, 2023). Notably, the impacts of those abiotic variable stressors vary by different seagrass species (Marbà et al., 1996; Lee et al., 2005; Celia and Fernando, 2023).

It is well known that atmospheric CO₂ levels are increasing under human influence, which pose many problems for marine ecosystems (Wang et al., 2022). The oceans continue to absorb CO₂ from the atmosphere and lead to a decline in ocean pH, changes in the carbonate system, and ocean acidification (OA) (Caldeira and Wickett, 2003), as a major disturbance in modern ocean events (Zhong et al., 2021). Environmental changes caused by OA will directly alter the function and structure of marine ecosystems, such as seagrasses (Gattuso et al., 2015; Zhong et al., 2021). By the end of the century, sea surface pCO₂ are projected to exceed 1000 μatm (IPCC, 2013). Currently, many studies have shown the variation of CO₂ concentration and temperature can produce a combined response in marine organisms (Nagelkerken and Munday, 2015; Watson et al., 2018; Brandenburg et al., 2019; Gao et al., 2019; Van Colen et al., 2020; Zhong et al., 2021). However, little is known about the adaptive responses of seagrasses to the OA and seasonal temperature fluctuations at different life stage. In recent decades, seagrass meadows due to human activities and climate change is a global concern, and their restoration has received much attention (Orth et al., 2006a). Seagrass seedling transplantation has been a useful tool for seagrass restoration (Pereda-Briones et al., 2018) due to its low damage on donor seagrass beds (Seddon, 2004; Orth et al., 2006b) and its ability to ensure the genetic diversity of seagrasses (Waycott et al., 2009). Therefore, it is important for seagrass conservation to explore the mechanisms by which seagrass seedlings respond to seasonal temperature increases and OA.

In the study, we studied the acclimatization response of Z. marina seedlings to OA and seasonal temperature changes. Under mesoscale experiments, this study investigated how Z. marina seedlings responded to (1) a large increase in seasonal temperature; (2) the occurrence of OA events; (3) OA events and seasonal temperature alternation, as assessed by growth, transcriptomics and metabolomics. This study aims to understand how seagrass seedlings adapt to seasonal temperature differences and global OA events in order to facilitate effective methods of seagrass conservation.

2 Materials and methods

2.1 Sample collection

A total of 120 Z. marina were collected in March 2022 at the sea area of Caofeidian (39.1131°N, 119.120°E), and 10 meters were spaced between each sample point to avoid collecting the same plants (Watanabe et al., 2005). All plants were in dark plastic boxes contaning in situ seawater to the laboratory within 1 h. The leaves, roots and stems of each plants were gently brushed and rinsed with tap water. Then, put them in a seawater recirculation tank for 7 days as a pre-culture stage. The transient condition was achieved in a recirculating seagrass culture system with a temperature of 14°C, a salinity of 32, a water flow rate of 800 mL min^-1, a light intensity of 240 µmol photons m^-2, and a 12:12 h light/dark ratio.

2.2 Experimental design

To setup temperature of each treatment, we used the data from the National Oceanic and Atmospheric Administration optimized interpolated sea surface temperature (OISST) (https://www.ncei.noaa.gov/metadata/geoportal/rest/metadata/item/gov.noaa.ncdc:C00844/html) (Richard et al., 2008), which spatial grid resolution is 0.25 degrees and temporal resolution is 1 day. The SST daily climatology of the Bohai Yellow Sea (34°N–42°N, 115°E–127°E) was extracted and calculated using the 20-year (2003–2022) SST mean values. All data were calculated in Ferret Version 7.6.0 and plotted in MATLAB software.

Four groups were set up with 30 well-grown, similarly growing plants in each group. Each group had a separate recirculation system consisting of a tank, a light system, a circulation system, a filtration system, a temperature control system and a monitoring system (more details in the Methods of Yan et al., 2023). Meanwhile, in situ seawater was used to maximize the simulation of the natural environment for Z. marina growth. The average SST of May is 14.03°C and the avarage SST of Jun is 19.02°C (Figure 1). To predict Z. marina response to OA, we choose the pCO₂ of the end of the century (1000 μatm) was used as the future OA condition. The control group (CK) was with the seawater temperature in May (14 ± 1°C) represented the spring season (Figure 1A), and the CO₂ concentration was that of the current atmospheric pCO₂ (410 μatm); the second group (OA) was the acidification group, in which the seawater temperature was 14 ± 1°C, and the CO₂ concentration was that of the future atmospheric pCO₂ (1000 ± 50 μatm); the third group (OW) was the SST temperature with seasonal changes treatment group, which has a current CO₂ concentration of 410 μatm and June temperature of 19 ± 1°C (Figure 1B); and the fourth group (AW) was the treatment group with composite stressors, i.e. extreme CO₂ concentration of 1000 ± 50 μatm and early summer temperature of 19 ± 1°C. During the experiment, the CO₂ flux of each group was controlled and monitored by a CO₂ enricher (SFO-E, China), which can pass CO₂ with a set concentration or air directly into the seawater tank (Yan et al., 2023). Water temperature, pH, and TDS were monitored in real-time using a temperature monitor (Caperplus Q2-PH/TDS, China) and a multiparameter water quality tester (ProQuatro; YSI, Yellow Springs, Ohio, USA) to make sure constant water quality. The whole treatment was 14 days.

FIGURE 1

Figure 1 Sea surface temperature daily climatology of May and June. (A) daily SST of May; (B) daily SST of June. The dashed black line denotes the average SST (blue line).

2.3 Plant growth

Nine samples were selected from each group. Since the plants did not have a very long leaf length (10 cm on average), a uniform de-averaging of the leaf length was not done. Leaf length (cm) was measured from the base of the meristematic tissue, leaf width (cm) was from the middle part of the leaf, and wet weight (g) was of the whole plant at the day 0 and 14 of the experiment. Before measurement, tap water was used to gently washed each sample, and then three growth index were measured to analyze their changes.

2.4 RNA extraction and transcriptome analysis

Three samples per treatment group were selected for transcriptomics sequencing. First, total RNA was extracted using the Trizol kit (Invitrogen, Carlsbad, CA, USA) according to the method provided by the manufacturer. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and assayed using an RNase-free RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies Palo Alto, CA, USA) and detected using RNase-free agarose gel electrophoresis. After extraction of total RNA, mRNA was enriched with Oligo(dT) beads, purified, and amplified by polymerase chain reaction. The resulting cDNA library was sequenced at Illumina Novaseq6000 to obtain all the raw data. To obtain high-quality clean data, the reads were first filtered using fastp V. 0.18.0 (Chen et al., 2019). Next, the reads were mapped using Bowtie2 V. 2.2.8 (Langmead and Salzberg, 2012), followed by removing the rRNA-mapped reads. Then, the remaining clean data were mapped to the reference genome of Z. marina using HISAT2. The reference gene from ENSEMBL (http://www.ensembl.org/index.html) was sequenced in 2016 and renewed their genome annotation in 2017 (Olsen et al., 2016). Fragments per kilobase of transcript per million mapped reads (FPKM) values were calculated using the RSEM software (Li and Dewey, 2011) to compare gene expression differences among samples. Principal component analysis (PCA) was performed using the R package gmodels (http://www.rproject.org/) to reveal relationships between samples. Differential expression analysis (Love et al., 2014) was performed on RNAs between two different groups using DESeq2 software (Robinson et al., 2010). Genes with a false discovery rate (FDR) of less than 0.05 and an absolute fold change of ≥ 2 were considered differentially expressed genes (DEGs). All DEGs were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/) to identify the major biological functions performed by DEGs. Kyoto Encyclopedia of the Genome (KEGG) metabolic pathway enrichment analysis was performed to identify major metabolic pathways in which DEGs were involved (Kanehisa and Goto, 2000).

2.5 Metabolite extraction and metabolome analysis

Six samples from each group were selected for metabolomics analysis. Aliquots of samples were prepared as quality control (QC) (Want et al., 2013). For metabolite extraction, we first ground the leaves (100 mg) individually in liquid nitrogen and then resuspended the homogenate in 500 µl of precooled 80% methanol by vortexing thoroughly. The samples were incubated on ice for 5 min and then centrifuged at 15,000 g for 20 min at 4°C. The supernatant was analyzed by LC-MS grade. The supernatant was diluted with LC-MS grade water to a final concentration of 53% methanol. The sample was then transferred to a new Eppendorf tube and centrifuged at 15,000 g for 20 minutes at 4°C. The supernatant was injected into the LC-MS sample. Finally, the supernatant was injected into the LC-MS/MS system for analysis (Dunn et al., 2011).

The metabolites were detected in positive ion mode (POS) and negative ion mode (NEG). Positive and negative ion modes were subsequently analyzed separately. Data were z-score normalized (Kolde, 2015). Orthogonal projection to latent structures-discriminant analysis (OPLS-DA) was performed on six comparison groups (CK-vs-AW, CK-vs-OW, OA-vs-OW, CK-vs-OA, OW-vs-AW, OA-vs-AW) using R package models (http://www.r-project.org/) (Bylesjö et al., 2006). The OPLS-DA model was further validated by cross-validation and permutation tests (Saccenti et al., 2014). The variable importance (VIP) score was applied to rank the most distinguishable metabolites between the two groups. The p values were calculated using R (Rao et al., 2016), and correlation heatmaps were plotted by the R corrplot software package (Wei et al., 2017). In addition, a t-test was used to screen for differentially expressed metabolites (p value < 0.05 and VIP ≥ 1). The abundance of differentially expressed metabolites in the same group was normalized by z-value using R, and the VIP scores from OPLS-DA were used to map the volcano plots (Yoon et al., 2020). Finally, metabolites were mapped to KEGG metabolic pathways enrichment analyses (Kanehisa and Goto, 2000). The p values were corrected for false discovery rate (FDR), and metabolic pathways with an FDR of ≤ 0.05 were screened as metabolic pathways significantly enriched for differentially expressed metabolites.

2.6 Statistical analysis

The data of growth indexes of each treatment were expressed as mean ± standard deviation (SD) and statistically analyzed. One-way analysis of variance was performed using IBM SPSS Statistics 27.0 (IBM SPSS, Chicago, IL, USA). Data were visualized using Origin 2021 (Originlab, Northampton, MA, USA).

3 Results

3.1 Plant growth

After 14 days, CK and OA tended to increase in wet weight at 14°C, and the increase was higher in the OA group; at 19°C, there was no significant change in OW and AW (Figure 2A). In addition, leaf length increased in all four groups (p < 0.05, Figure 2B). The increase was smaller in OA, OW, and AW groups than in CK group. OW group had the longest leaf length, which showed the largest increase, while the AW group showed the smallest. Similarly, leaf width increased significantly in all four groups (Figure 2C, p < 0.05), with the largest increase in CK group, followed by OA group, and the samllest increase in AW group.

FIGURE 2

Figure 2 (A) Wet weight, (B) leaf length, and (C) leaf width for different groups at the beginning and end of the experiment. CK means the group control with 14°C, and current atmospheric pCO₂ (410 μatm); OA is with seawater temperature of 14°C, and the future atmospheric pCO₂ (1000 μatm); OW is with a current CO₂ concentration (410 μatm) and June temperature of 19°C; AW is with an extreme CO₂ concentration of 1000 μatm and early summer temperature of 19°C.

3.2 Transcriptome results

A total of 12 leaf tissue samples (three replicates from each group) were sequenced and got a 9.55 Gb high-quality reads with an average of 4.96 Gb. The percentage of bases filtered to Q30 in clean data ranged from 92.56% to 94.53%. The GC ratios of filtered sequence bases ranged from 44.79% to 45.63%. Comparison of the clean data with the Z. marina reference genome (Olsen et al., 2016) revealed a rate of 90.25–94.57%.

Based on transcriptome analysis, we found that the number of DEGs varied among different comparisons (Figure S1). The greatest number of DEGs were found between CK and AW. This was followed by CK-VS-OW, OA-VS-OW, CK-VS-OA, OW-VS-AW, and the least was OA-VS-AW.

A total of 1,290 DEGs were identified in AW compared to CK, of which 763 were upregulated genes (Figure 3A). GO enrichment analysis showed that DEGs were mainly enriched in biological process, cellular process, metabolic process, and so on. The top 10 terms in biological process were flavonoid biosynthetic process, asparagine biosynthetic process, asparagine metabolic process, flavonoid metabolic process, glutamine family amino acid metabolic process, single-organism process, hyperosmotic response, response to chemical (q < 0.01). These DEGs were involved in 107 KEGG pathways (Figure S2A). Among them, biosynthesis of secondary metabolites; flavonoid biosynthesis; alanine, aspartate and glutamate metabolism; glycine, serine and threonine metabolism; biosynthesis of amino acids, starch and sucrose metabolism; and galactose metabolism showed significant enrichment (p < 0.01, Figure S2A), with the following valine, leucine and isoleucine biosynthesis; arachidonic acid metabolism; DNA replication; plant-pathogen interaction; betalain biosynthesis; circadian rhythm-plant, 2-oxocarboxylic acid metabolism; and vitamin B6 metabolism (p < 0.05, Figure S2A).

FIGURE 3

Figure 3 Volcano plot of differentially expressed genes (DEGs) in different treatments. DEGs in (A) CK-vs-AW, (B) CK-vs-OW, (C) OA-vs-OW, (D) CK-vs-OA, (E) OW-vs-AW, (F) OA-vs-AW. Each dot is a gene. The horizontal coordinate is log₂ fold-change, and the vertical coordinates is -log10 q values. Red, significantly upregualted genes; blue, significantly downregualted genes; gray, no significant differences. All DEGs are with q < 0.05 and |log₂ FC| > 1. The closer the genes are to the ends, the greater the degree of difference. CK means the group control with 14°C, and current atmospheric pCO₂ (410 μatm); OA is with seawater temperature of 14°C, and the future atmospheric pCO₂ (1000 μatm); OW is with a current CO₂ concentration (410 μatm) and June temperature of 19°C; AW is with an extreme CO₂ concentration of 1000 μatm and early summer temperature of 19°C.

In OW compared to CK, 1,220 DEGs were identified with 701 upregulated genes (Figure 3B). GO enrichment analysis showed that the DEGs significantly enriched in biological processes were response to stress, chemical, organic substance, cellular response to UV (q < 0.005), response to stimulus and acid chemical (q < 0.05). KEGG pathways revealed these DEGs were enriched in 103 pathways (Fugire S2b). The most significant pathways were metabolic pathways, amino sugar and nucleotide sugar metabolism, starch and sucrose metabolism, plant-pathogen interaction (p < 0.01, Figure S2B), followed by biosynthesis of secondary metabolites, arachidonic acid metabolism, photosynthesis-antenna proteins, galactose metabolism, nitrogen metabolism, flavonoid biosynthesis, DNA replication, fructose and mannose metabolism, ascorbate and aldarate metabolism, valine, leucine and isoleucine biosynthesis (p < 0.05, Figure S2B). In particular, NCED (9-cis-epoxycarotenoid dioxygen) was significantly upregulated in carotenoid biosynthesis.

A total of 1146 DEGs were identified in OW compared to OA, of which 555 were upregulated and 591 were downregulated (Figure 3C). GO enrichment analysis showed that response to organic substance, stimulus, acid chemical, chemical, endogenous stimulus, stress, generation of precursor metabolites and energy, response to organonitrogen compound, photosynthesis, light reaction, photosynthesis, response to hormone, asparagine biosynthetic process, cellular response to acid chemical, and immune system process were significantly enriched in biological processes (q < 0.005). There were 92 significant KEGG pathways (Figure S2C). Among them, metabolic pathways were the most significant (p < 0.01 and q < 0.005), followed by biosynthesis of secondary metabolites; photosynthesis-antenna proteins; glutathione metabolism; carbon fixation in photosynthetic organisms; pentose phosphate pathway; photosynthesis; MAPK signaling pathway-plants; alanine, aspartate and glutamate metabolism; and biosynthesis of amino acid (p < 0.01).

A total of 836 DEGs were identified in OA compared to CK, of which 442 were upregulated (Figure 3D). In biological processes, DEGs significantly enriched in, firstly, asparagine biosynthetic process, asparagine metabolic process, response to organonitrogen compound, aspartate family amino acid biosynthetic process, response to chemical, glutamine family amino acid metabolic process, response to stimulus, and acid chemical (q < 0.001), followed by response to stress, cell communication, response to light intensity, cellular response to starvation, response to starvation, cellular response to nutrient levels, hyperosmotic response, alpha-amino acid biosynthetic process, response to nutrient levels, cellular response to extracellular stimulus, cellular response to external stimulus, response to extracellular stimulus (q < 0.005). These DEGs were involve in 98 KEGG pathways (Figure S2D). Specifically speaking, biosynthesis of secondary metabolites; alanine, aspartate and glutamate metabolism; biosynthesis of amino acids; carbon metabolism; carbon fixation in photosynthetic organisms; arachidonic acid metabolism; glycine, serine and threonine metabolism; linoleic acid metabolism; and phenylalanine, tyrosine and tryptophan biosynthesis showed the significance (p < 0.01, Figure S2D), followed by starch and sucrose metabolism, glyoxylate and dicarboxylate metabolism, TCA cycle, cyanoamino acid metabolism, vitamin B6 metabolism, flavonoid biosynthesis, alpha-linolenic acid metabolism, and galactose metabolism (p < 0.05, Figure S2D). In tropane, piperidine and pyridine alkaloid biosynthesis, hisC (histidinol-phosphate aminotransferase) was significantly upregulated. In starch and sucrose metabolism, eight genes showed strong up or downregulation trends. For example, pgm (phosphoglucomutase) and glgA (starch synthase) showed upregulation, and INV, sacA (beta-fructofuranosidase), and SUS (sucrose synthase) were downregulated.

A total of 591 DEGs were identified in AW compared to OW, of which 380 were upregulated genes (Figure 3E).Top 15 GO terms in biological process were response to stimulus, chemical, organic substance, acid chemical, endogenous stimulus, organonitrogen compound, signal transduction, hormone, other organisms, biotic stimulus, external biotic stimulus, stress, cellular response to stimulus, signaling, and single-organism signaling (q < 0.001). KEGG results showed that DEGs were enriched in 78 metabolic pathways. Biosynthesis of secondary metabolites; amino sugar and nucleotide sugar metabolism; phenylpropanoid biosynthesis; flavonoid biosynthesis; phenylalanine, tyrosine and tryptophan biosynthesis; vitamin B6 metabolism; MAPK signaling pathway-plant, ascorbate and aldarate metabolism; glycine, serine and threonine metabolism; betalain biosynthesis; and tyrosine metabolism expressed significantly (p < 0.01, Figure S2E). Isoquinoline alkaloid biosynthesis; starch and sucrose metabolism; cutin, suberine and wax biosynthesis; linoleic acid metabolism; and biosynthesis of amino acids likewise showed a significance (p < 0.05, Figure S2E).

There were 172 DEGs in AW compared to OA, of which 132 were upregulated (Figure 3F). Protein complex subunit organization and protein oligomerization were significantly enriched in biological processes (p < 0.05). DEGs were enriched in 47 metabolic pathways. DNA replication, photosynthesis-antenna proteins, RNA degradation, sesquiterpenoid and triterpenoid biosynthesis, and carbon fixation in photosynthetic organisms showed significant enrichment (p < 0.05, Figure S2F).

3.3 Metabolomic results

We selected 24 samples (6 samples per group) for metabolomics analysis. Two ion modes, positive ion mode (POS) and negative ion mode (NEG), were combined, and all data were analyzed separately in each mode. A total of 25,871 metabolites were identified and 6,028 were named metabolites. Under POS mode, 3,723 named metabolites were identified, mainly including glycerolipids (15%), carboxylic acids and derivatives (9%), and organooxygen compounds (8%). A total of 2,305 named metabolites composed of glycerolipids (17%), glycerophospholipids (9%), and fatty acyls (8%)were identified under NEG mode. There were differences among different comparisons (Figure 4). The prediction rates of the OPLS-DA models were 0.413-0.794 (Table S2), which means models were reliable (Figure S3).

FIGURE 4

Figure 4 Score plot of orthogonal partial least-squares discriminant analysis (OPLS-DA). (A-F): POS mode; (G-L): NEG mode. (A, G): CK-vs-AW; (B, H): CK-vs-OW; (C, I): OA-vs-OW; (D, J): CK-vs-OA; (E, K): OW-vs-AW; (F, L): OA-vs-AW. The horizontal coordinate represents the predicted principal component (PC), which can be used to directly distinguish differences between different treatments; the vertical coordinate represents the orthogonal PC, which can reflect the variation within one group. CK means the group control with 14°C, and current atmospheric pCO₂ (410 μatm); OA is with seawater temperature of 14°C, and the future atmospheric pCO₂ (1000 μatm); OW is with a current CO₂ concentration (410 μatm) and June temperature of 19°C; AW is with an extreme CO₂ concentration of 1000 μatm and early summer temperature of 19°C.

There were differences in the differentially expressed metabolites between different comparisons (Figure S4). In AW compared with CK, 114 differentially expressed metabolites were identified in POS mode, 69 were upregulated, and the top five were kojibiose, D-maltose, 4-hydroxy-2-butenoic acid gamma-lactone, 1,2,3-trihydroxybenzene, and 3-hydroxyadipic acid 3,6-lactone (Figure 5). In NEG mode, 35 differentially expressed metabolites were identified, of which 28 were upregulated (Figure S4). Among them, 1-deoxy-a-(N6-lysino)-D-fructose, 2-Ethyl-2-hydroxybutyric acid, malonic acid, L-glutamic acid, and itaconic acid were the top five, and only N-acetylglutamic acid was downregulated (Figure 5). KEGG results revealed 67 pathways were significantly enriched (Figure 6A), of which biosynthesis of plant secondary metabolites and biosynthesis of terpenoids and steroids were the most highly enriched (p < 0.01), followed by the butanoate metabolism; catecholamine transferase inhibitors; propanoate metabolism; biosynthesis of alkaloids derived from ornithine, lysine and nicotinic acid; C5-Branched dibasic acid metabolism; sulfur metabolism; and biosynthesis of alkaloids derived from shikimate pathway (p < 0.05).

FIGURE 5

Figure 5 Variable importance in projection (VIP) plots of differentially expressed metabolites in different treatments. (A-F): POS mode; (G-L): NEG mode. (A, G): CK-vs-AW; (B, H): CK-vs-OW; (C, I): OA-vs-OW; (D, J): CK-vs-OA; (E, K): OW-vs-AW; (F, L): OA-vs-AW. The horizontal coordinate is the VIP value, and the vertical coordinate is the top 15 differentially expressed metabolites; the abundance of each metabolite in each group was averaged, and then z-score analysis was performed for each group. The colors on the right side indicated the abundance of the metabolite in different groups, red indicated upregulation, and green indicated downregulation. The higher the value of VIP, the greater their contribution to the differentiation of samples. The metabolites with default VIP greater than 1 had a significant difference. CK means the group control with 14°C, and current atmospheric pCO₂ (410 μatm); OA is with seawater temperature of 14°C, and the future atmospheric pCO₂ (1000 μatm); OW is with a current CO₂ concentration (410 μatm) and June temperature of 19°C; AW is with an extreme CO₂ concentration of 1000 μatm and early summer temperature of 19°C.

FIGURE 6

Figure 6 Bubble diagram of metabolic pathways. Top 20 pathways with the smallest q values were shown. The vertical coordinate is the pathway name, the horizontal coordinate is the enrichment factor, and the cicle size indicates how many metabolites in the pathway. The color shows the significance, and the redder the color, the smaller the q value. CK means the group control with 14°C, and current atmospheric pCO₂ (410 μatm); OA is with seawater temperature of 14°C, and the future atmospheric pCO₂ (1000 μatm); OW is with a current CO₂ concentration (410 μatm) and June temperature of 19°C; AW is with an extreme CO₂ concentration of 1000 μatm and early summer temperature of 19°C.

A total of 109 differentially expressed metabolites were identified in OW compared with CK under POS mode, including 71 upregulated. Among them, L-glutamine (downregulated), kojibiose, 4-hydroxy-2-butenoic acid gamma-lactone, 1,2,3-Trihydroxybenzene, and D-maltose were the top five (Figure 5). In NEG mode, there were 27 differentially expressed metabolites with 22 upregulated and five downregulated (Figure S4). And cis,cis-Muconic acid, trehalose, 2-ketobutyric acid, succinic acid, malonic acid, and D-ribose were differentially expressed metabolites, and all were upregulated (Figure 5). KEGG pathway enrichment analysis revealed 65 noteworthy signaling pathways (Figure 6B), among which plant hormone signal transduction carotenoid biosynthesis was most highly enriched (p < 0.01), followed by aminoacyl-tRNA biosynthesis; alanine, aspartate and glutamate metabolism; catecholamine transferase inhibitors; biosynthesis of plant secondary metabolites; and biosynthesis of terpenoids and steroids (p < 0.05). Beta-cryptoxanthin, zeaxanthin, and abscisate were significantly upregulated in carotenoid biosynthesis.

In OW compared with OA, 131 differentially expressed metabolites were identified in POS mode, of which 82 were upregulated (Figure S4), and L-glutamine, L-isoleucine, 6-deoxyfagomine, and kojibiose (downregualted) were the top five (Figure 5). In NEG mode, 48 were identified, of which 29 were upregulated (Figure S4), and cis,cis-Muconic acid; 3,4-Dihydroxybenzaldehyde;1-deoxy-1-(N6-lysino)-D-fructose (only downregulated); succinic acid; 2-ketobutyric acid; and caffeic acid were top five (Figure 5). KEGG analysis revealed 86 significant pathways (Figure 6C), of which aminoacyl-tRNA biosynthesis; valine, leucine, and isoleucine biosynthesis; biosynthesis of plant secondary metabolites; glucosinolate biosynthesis; ABC transporter (p < 0.01); propanoate metabolism; biosynthesis of amino acid; glycine, serine and threonine metabolism; and catecholamine transferase inhibitors (p < 0.05) were significantly enriched.

In OA compared with CK, 65 differentially expressed metabolites were identified in POS mode, of which 37 were upregulated (Figure S4). The top five were L-glutamine, 6-deoxyfagomine, beta-leucine, betaine, and piperidine and all of them showed upregulation (Figure 5). A total of 21 were identified in NEG mode, of which 12 were upregulated (Figure S4), and 1-deoxy-1-(N6-lysino)-D-fructose, adenine, glyceraldehyde (downregulated), D-ribose, succinic acid (downregulated) were the highly significant differentially expressed metabolites (Figure 5). KEGG pathway analysis of differentially expressed metabolites identified 61 significant pathways (Figure 6D), such as glucosinolate biosynthesis; plant hormone signal transduction; aminoacyl-tRNA biosynthesis; catecholamine transferase inhibitors (p < 0.01); biosynthesis of plant secondary metabolites; ABC transporters; starch and sucrose metabolism; tropane, piperidine and pyridine alkaloid biosynthesis (p < 0.05). Among them, L-isoleucine and azacyclohexane were significantly upregulated in tropane, piperidine, and pyridine alkaloid biosynthesis. In starch and sucrose metabolism, maltose showed upregulation, and trehalose 6-phosphate showed downregulation.

In AW compared with OW, 63 differentially expressed metabolites were identified in POS mode, of which 29 were upregulated metabolites; and 26 were identified in NEG mode, of which 14 were upregulated metabolites (Figure S4). In POS mode, L-glutamine, 6-deoxyfagomine, N-acryloylglycine, piperidine, and betaine were the top five, and all of them showed upregulation (Figure 5). In NEG mode, cis,cis-muconic acid (downregulated), 1-deoxy-1-(N6-lysino)-D-fructose, succinic acid (downregulated), and oleamide were highly significant (Figure 5). There were 41 significant KEGG pathways (Figure 6E), such as alanine, aspartate and glutamate metabolism; sulfur metabolism (p < 0.01); catecholamine transferase inhibitors, oxidative phosphorylation; tyrosine metabolism; butanoate metabolism; propanoate metabolism; degradation of aromatic compounds; glyoxylate and dicarboxylate metabolism; lysine degradation; nitrogen metabolism; and phenylalanine metabolism (p < 0.05). L-glutamine was significantly upregulated in nitrogen metabolism.

In AW compared with OA, 64 differentially expressed metabolites were identified in POS mode, of which 36 were upregulated, and 26 were identified in NEG mode, of which 14 were upregulated, (Figure S4). In POS mode, kojibiose, D-maltose, and 3-hydroxyadipic acid 3,6-lactone were the top three, and all of them showed upregulation (Figure 5). In NEG mode, succinic acid, eicosapentaenoic acid (downregulated), and oleamide were highly significant differentially expressed metabolites (Figure 5). KEGG results identified 59 noteworthy signaling pathways (Figure 6F), among which biosynthesis of alkaloids derived from histidine and purine, type I polyketide structures, butanoate metabolism, propanoate metabolism, lysine degradation, phenylalanine metabolism, biosynthesis of plant hormones, glyoxylate and dicarboxylate metabolism, pyruvate metabolism, TCA cycle, sulfur metabolism, carbon metabolism, and biosynthesis of plant secondary metabolites were the most enriched (p < 0.01), followed by biosynthesis of alkaloids derived from terpenoid and polyketide; fatty acid degradation; oxidative phosphorylation; biosynthesis of alkaloids derived from ornithine, lysine and nicotinic acid; biosynthesis of terpenoids and steroids; biosynthesis of alkaloids derived from shikimate pathway; GABA-A receptor agonists/antagonists; valine, leucine and isoleucine degradation; and degradation of aromatic compounds (p < 0.05). Significant upregulation of acetyl-CoA and succinate occurred in TCA cycle.

4 Discussion

It is well known that seagrass beds have important conservation values (Barbier et al., 2011). However, the response of Z. marina to the single and synergistic effects of seasonal warming and OA, especially at the seedling stage, is limited. In this study, we explored the mechanisms and response differences in different environmental combinations (OW, OA, AW) to Z. marina seedling growth (leaf length, leaf width, wet weight), RNAs and metabolites from multiple perspectives using phenomics, transcriptomics and metabolomics.

4.1 Seasonal temperature increase could limit the growth of Zostera marina seedlings

The results of metabolomics and transcriptomics showed that under suitable temperature, Z. marina seedlings activate the upregualtion of beta-cryptoxanthi, zeaxanthin, and abscisate upregualted by NCED and then ultimately promotes the biosynthesis of carotenoid biosynthesis to adapt to rapid seasonal temperature increases (Figure 7). Notably, this pathway was specific to such a condition (Figure 7). Previous studies have shown that several carotenoids at different depths exhibit stimulatory properties in Posidonia sinuos (Collier et al., 2008). The process of epoxidation of zeaxanthin requires oxygen and energy consumption (Janhs et al., 2009), which may lead to reduced photosynthetic efficiency and competition for energy with other metabolic processes (Silva et al., 2013). This is detrimental to Z. marina seedlings under seasonally rapid SST warming. Also, the results of the multiple-omics may correspond to the results of growth. After 14 days, there was no significant increase in wet weight in OW group. Although the leaf length and width increased, the growth in OW was not as good as that of CK (Figure 3). This also suggests that, for seedlings in group OW, although seasonal temperature increases over a short period are not fatal to individual plants, they can still limit the growth of Z. marina seedlings for a certain period. Although many studies have shown that summer is the peak growth season for Z. marina (Nguyen et al., 2021; Azcarate-Garcia et al., 2022), this is different from results of this study. This may be due to the fact that our study was the result of intense stress for a short period. Since temperature increase occurred over a longer period, Z. marina seedlings gradually adapted to the change in temperature. In addition, as Z. marina grow from seedlings to adults, their ability to adapt to temperature increases. At this point, the temperature of 19°C is no longer a stressor that inhibits the growth of Z. marina adult (Nguyen et al., 2021; Azcarate-Garcia et al., 2022; Yan et al., 2023).

FIGURE 7

Figure 7 Kyoto encyclopedia of genes and genomes pathways based on p-value. Different classes are shown in different colors. The redder the color, the smaller the p-value. CK means the group control with 14°C, and current atmospheric pCO₂ (410 μatm); OA is with seawater temperature of 14°C, and the future atmospheric pCO₂ (1000 μatm); OW is with a current CO₂ concentration (410 μatm) and June temperature of 19°C; AW is with an extreme CO₂ concentration of 1000 μatm and early summer temperature of 19°C.

4.2 Ocean acidification may promote the growth of Zostera marina seedlings under suitable temperature

Under OA, by upregulating L-isoleucine and azacyclohexane, Z. marina seedlings promote the biosynthesis of tropane, piperidine, and pyridine alkaloids (TPPA). Meanwhile, starch and sucrose metabolism were promoted through the coordination of several genes, such as upregualtion of glgA and downregualtion of sucrose synthase, and through the regulation of some metabolites, such as the upregualtion of maltose and downregualtion of trehalose 6-phosphate (Figure 7). Finally, a different response mechanism from the other stress environments was realized (the two pathways mentioned above were not significant in the other groups) (Figure 7). During the growth of Z. marina seedlings, leaf length, leaf width and wet weight increased significantly, with the highest increase in wet weight compared to the other comparisons. This suggests that OA did not inhibit the growth of Z. marina seedlings at suitable temperatures and even played a role in promoting it (Invers et al., 2002; Jiang et al., 2010; Russell et al., 2013; Martínez-Crego et al., 2014; Campbell and Fourqurean, 2018; Guerrero-Meseguer et al., 2020). In many higher plants, low molecular weight chemicals such as alkaloids are produced, and these secondary metabolites provide immunity to the body against pathogens (Yamada and Sato, 2013). The TPPA biosynthetic pathway, which can be used to modulate alkaloids with pharmacological properties such as anti-inflammatory and antimicrobial properties, has been frequently reported in previous studies (Rischer et al., 2006; Huang et al., 2019). We hypothesized that biosynthesis of TPPA, as well as starch and sucrose metabolism, play distinct roles in metabolic capacity from the other comparisons (Figure 7), thereby increasing the resistance of Z. marina seedlings and enhancing their adaptation to OA.

In summary, under suitable temperatures, both in the face of increased temperature and OA, biosynthesis of plant secondary metabolites, plant hormone signal transduction, catecholamine transferase inhibitors, aminoacyl-tRNA biosynthesis, all of which enhanced the biosynthesis of plant secondary metabolites in Z. marina seedlings (Figure 7). Among them, plant hormone signal transduction and aminoacyl-tRNA biosynthesis were unique when Z. marina seedlings responded to temperature increases or OA (Figure 7).

4.3 Sequential or concurrent occurence of seasonal temperature increase and ocean acidification differentially affect Zostera marina seedlings

In OW-VS-AW, nitrogen metabolism showed unique significance compared with other metabolic pathways (Figure 7). That is, under high SST caused by seasonal change, Z. marina seedlings increase the ability to adapt to OA by upregulating L-glutamine to accelerate nitrogen metabolism in such a unique metabolic pathway. In OA-VS-AW, we also found significant enrichment of metabolic pathways different from the other treatments, including valine, leucine and isoleucine degradation; TCA; pyruvate metabolism; biosynthesis of alkaloids derived from histidine and purine; biosynthesis of alkaloids derived from terpenoid and polyketide; and biosynthesis of plant hormones (Figure 7). Notably, TCA cycle can link almost individual metabolic pathways (Akram, 2014; Lin et al., 2016; Zhang et al., 2017), which is an important central pathway can energize organisms and produce essential metabolites for energy production and biosynthetic reaction (Arnold et al., 2022). TCA cycle did not show a significant effect in CK-VS-OW, although they both increased the temperature (Figure 7). This suggests that the stress mechanisms of Z. marina seedlings are different even if they later encounter the same environmental stress. Under OA, Z. marina seedlings respond to seasonal temperature surges by a unique metabolic pathway that promotes TCA cycle through the upregualtion of acetyl-CoA and succinate (Arnold et al., 2022). In addition to these treatment-specific KEGG pathways, there are co-existing metabolic pathways that counteract high temperatures and OA, such as lysine degradation, phenylalanine metabolism, glyoxylate and dicarboxylate metabolism, oxidative phosphorylation, and degradation of aromatic compounds.

These results suggested Z. marina seedlings respond to environmental changes by accelerating common metabolic pathways when faced with SST increase by seasonal change and OA. However, it is worth noting that the sequencing of these two events had different effects on Z. marina seedlings. In addition, OW-VS-AW and CK-VS-OA, and OA-VS-AW and CK-VS-OW also shared the same metabolic pathways but with many differences, which means although Z. marina responds to the same oceanic events, their initial growth environment has some influence on their response mechanisms.

Metabolomics results showed that Z. marina seedlings responds to the dual stress conditions of seasonal SST increases and OA simultaneously by actively mobilizing nine KEGG pathways, such as the biosynthesis of plant secondary metabolites, biosynthesis of terpenoids and steroids, carbon metabolism, butanoate metabolism (Figure 7). Among them, propionate metabolism, butyrate metabolism, and sulfur metabolism showed enrichment under different stress conditions (CK-VS-AW, OA-VS-AW, and OW-VS-AW) (Figure 7). This suggests that, regardless of whether the increases in temperature and OA occurred sequentially or concurrently, Z. marina increased these three metabolisms in order to improve their adaptive capacity to environmental changes. In terms of growth of Z. marina seedlings, AW had the lowest growth rate compared with other treatments, and the wet weight did not increase basically. Compared with single stress, OA had the least effect on the growth, followed by OW; AW had the most inhibitory effect. This result is inconsistent with the results of Yan et al., who believed that the addition of OA would reduce the effect of temperature increase on the growth of Z. marina (Yan et al., 2023). We think the reason for this difference is that the materials selected for the two studies had different growth periods of Z. marina. We selected the period of seedling, during which the plants are more sensitive to changes in environmental factors (Gogna and Bhatla, 2020). Transcriptomics and metabolomics were combined to show that the sequence of environmental factors greatly affected the response mechanism and adaptive ability of seagrass plants.

5 Conclusion

In this study, we investigated the response of Z. marina seedlings to seasonal warming and OA using phenomics and multi-omics (transcriptomics and metabolomics). Different combination or occurance order of environmental changes prompted different responces for Z. marina seedlings. Leaf length, leaf width, and wet weight increased significantly under different stress, with overall plant growth reaching an optimum under OA, but OW and AW were lower than the control. This indicated that OA did not inhibit the growth of plant seedlings and even promoted their growth. The short-term seasonal temperature increase did not cause fatal damage but still limited their growth to some extent. Multi-omics results showed that at appropriate temperatures, Z. marina seedlings collectively promote carotenoid biosynthesis through the upregualtion of beta-cryptoxanthin, zeaxanthin, and the upregualtion of NCED, which in turn mobilizes the upregualtion of abscisate, ultimately differentiating the treatment from other treatments for resistance to the stress mechanisms against seasonal warming. Under the same conditions, the plants promoted tropane, piperidine and pyridine alkaloid biosynthesis through the upregualtion of L-isoleucine azacyclohexane. The different genes (upregualtion of glgA and downregualtion of sucrose synthase) and the differentially expressed metabolites (upregualtion of maltose and downregualtion of trehalose 6-phosphate) were used to promote starch and sucrose metabolism, and finally, the specific mechanism to respond to the stress factors was realized to enhance the adaptive capacity to OA. When Z. marina seedlings are already under some environmental stress, it responds differently to seasonal warming and OA. Under high temperatures, when Z. marina seedlings respond to OA, they accelerate nitrogen metabolism by upregulating L-glutamine. In OA, nine unique metabolic pathways, including TCA cycle and pyruvate metabolism, work together to respond to environmental changes when ambient temperatures rise again. In response to two stressors at the same time, plants mobilize biosynthesis of plant secondary metabolites, carbon metabolism, and other metabolic pathways to increase their resistance. Based on these results, we found that although plants face the same factors, their initial survival environment leads to different response mechanisms. In addition to the specific pathways described above for each treatment, there are common pathways that stimulate plant metabolism and sulfur metabolism in response to seasonal temperature rise or/and OA, such as propionate metabolism, butyrate metabolism, and sulfur metabolism.

Data availability statement

The transcriptome data can be found in Bioproject accession number: PRJNA1031782. Further inquiries can be directed to the corresponding authors.

Author contributions

ZW: Formal Analysis, Methodology, Writing – original draft, Writing – review & editing, Data curation, Visualization. YP: Formal Analysis, Methodology, Writing – original draft, Data curation. WY: Conceptualization, Supervision, Writing – review & editing, Writing – original draft. LL: Conceptualization, Methodology, Writing – original draft. BZ: Conceptualization, Funding acquisition, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. The study was supported by the National Natural Science Foundation of China (grant number: 42076147). The funder is the corresponding author BZ.

Acknowledgments

We thank editors and reviewers for very helpful comments on the manuscript. Thanks for the language-editing service provided by the Charlesworth Group.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2023.1304132/full#supplementary-material

References

Akram M. (2014). Citric acid cycle and role of its intermediates in metabolism. Cell Biochem. Biophys. 68 (3), 475–478. doi: 10.1007/s12013-013-9750-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Arnold P. K., Jackson B. T., Paras K. I., Brunner J. S., Hart M. L., Newsom O. J., et al. (2022). A non-canonical tricarboxylic acid cycle underlies cellular identity. Nature 603 (7901), 477–481. doi: 10.1038/s41586-022-04475-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Azcarate-Garcia T., Beca-Carretero P., Cara C. L., Villamayor B., Cosnett E., Bermejo R., et al. (2022). Seasonal plant development and meadow structure of Irish and southern Spanish seagrass populations. Aquat. Botany 183, 4. doi: 10.1016/j.aquabot.2022.103569

CrossRef Full Text | Google Scholar

Baden S. P., Pihl L. (1984). Abundance, biomass and production of mobile epibenthic fauna in Zostera marina (L.) meadows, western Sweden. Ophelia 23 (1), 65–90. doi: 10.1080/00785236.1984.10426605

CrossRef Full Text | Google Scholar

Barbier E. B., Hacker S. D., Kennedy C., Koch E. W., Stier A. C., Silliman B. R. (2011). The value of estuarine and coastal ecosystem services. Ecol. Monographs. 81 (2), 169–193. doi: 10.1890/10-1510.1

CrossRef Full Text | Google Scholar

Boström C., Baden S., Krause-Jensen D. (2003). “The seagrasses of Scandinavia and the Baltic sea,” in World Atlas of Seagrasses: Present Status and Future Conservation. Eds. Green E., Short F. T. (California: California Press), 27–37.

Google Scholar

Brandenburg K. M., Velthuis M., Waal D. B. V. D. (2019). Meta-analysis reveals enhanced growth of marine harmful algae from temperate regions with warming and elevated CO₂ levels. Glob Chang Biol. 25, 2607–2618. doi: 10.1111/gcb.14678

PubMed Abstract | CrossRef Full Text | Google Scholar

Burke M. K., Dennison W. C., Morre K. A. (1996). Non-structural carbohydrate reserves of eelgrass Zostera marina. Mar. Ecol. Prog. Series. 137, 195–201. doi: 10.3354/meps137195

CrossRef Full Text | Google Scholar

Bylesjö M., Rantalainen M., Cloarec O., Nicholson J. K., Elaine Holmes E., Johan Trygg J. (2006). OPLS discriminant analysis: combining the strengths of PLS-DA and SIMCA classification. J. Chemometrics: A J. Chemometrics Soc. 20 (8-10), 341–351. doi: 10.1002/cem.1006

CrossRef Full Text | Google Scholar

Caldeira K., Wickett M. K. (2003). Anthropogenic carbon and ocean pH. Nature 425, 365. doi: 10.1038/425365a

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell J. E., Fourqurean J. W. (2018). Does nutrient availability regulate seagrass response to elevated CO₂? Ecosyst. Springer US. 21, 1269–1282. doi: 10.1007/s10021-017-0212-2

CrossRef Full Text | Google Scholar

Celia D. S., Fernando T. (2023). Annual variation in leaf canopy across a seagrass distribution range is modulated by local environmental variation. Mar. Biol. Res. 19 (2-3), 132–140. doi: 10.1080/17451000.2023.2193898

CrossRef Full Text | Google Scholar

Coll M., Schmidt A., Romanuk T., Lotze H. K. (2011). Food-Web structure of seagrass communities across different spatial scales and human impacts. PLoS One 6 (7), e22591–e22513. doi: 10.1371/journal.pone.0022591

PubMed Abstract | CrossRef Full Text | Google Scholar

Collier C. J., Lavery P. S., Ralph P. J., Masini R. J. (2008). Physiological characteristics of the seagrass Posidonia sinuosa along a depth-related gradient of light availability. Mar. Ecol. Prog. Ser. 353, 65–79. doi: 10.3354/meps07171

CrossRef Full Text | Google Scholar

Duarte C. M., Middelburg J. J., Caraco N. (2005). Major role of marine vegetation on the oceanic carbon cycle. Biogeosciences (online). 2 (1), 1–8. doi: 10.5194/bg-2-1-2005

CrossRef Full Text | Google Scholar

Dunn W. B., Broadhurst D., Begley P., Zelena E., Francis-McIntyre S., Anderson N., et al. (2011). Human Serum Metabolome (HUSERMET) Consortium. Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry. Nat. Protoc. 6 (7), 1060–1083. doi: 10.1038/nprot.2011.335

PubMed Abstract | CrossRef Full Text | Google Scholar

Espino F., Tuya F., Brito A., Haroun R. J. (2011). Ichthyofauna associated with Cymodocea nodosa meadows in the Canarian Archipelago (central eastern Atlantic): Community structure and nursery role. Cienc. Marinas. 37 (2), 157–174. doi: 10.7773/cm.v37i2.1720

CrossRef Full Text | Google Scholar

Fourqurean J. W., Duarte C. M., Kennedy H., Marbà N., Holmer M., Mateo M. A., et al. (2012). Seagrass ecosystems as a globally significant carbon stock. Nat. Geoscience. 5 (7), 505–509. doi: 10.1038/ngeo1477

CrossRef Full Text | Google Scholar

Gao K., Beardall J., Häder D. P., Hall-Spencer J. M., Gao G., Hutchins D. A. (2019). Effects of ocean acidification on marine photosynthetic organisms under the concurrent influences of warming, UV radiation, and deoxygenation. Front. Mar. Sci. 6, 322. doi: 10.3389/fmars.2019.00322

CrossRef Full Text | Google Scholar

Gattuso J. P., Magnan A., Billé R., Cheung W., Howes E., Joos F., et al. (2015). Contrasting futures for ocean and society from different anthropogenic CO₂ emissions scenarios. Science 349, aac4722. doi: 10.1126/science.aac4722

PubMed Abstract | CrossRef Full Text | Google Scholar

Gogna M., Bhatla S. C. (2020). Salt-tolerant and -sensitive seedlings exhibit noteworthy differences in lipolytic events in response to salt stress. Plant Signal Behav. 15 (4), 1737451. doi: 10.1080/15592324.2020.1737451

PubMed Abstract | CrossRef Full Text | Google Scholar

Granata T. C., Serra T., Colomer J., Casamitjana X., Duarte C. M., Gacia E. (2001). Flow and particle distributions in a nearshore seagrass meadow before and after a storm. Mar. Ecol. Prog. Series. 218, 95–106. doi: 10.3354/meps218095

CrossRef Full Text | Google Scholar

Guerrero-Meseguer L., Cox T. E., Sanz-Lázaro C., Schmid S., Enzor L. A., Major K., et al. (2020). Does ocean acidification benefit seagrasses in a mesohaline environment? A mesocosm experiment in the Northern Gulf of Mexico. Estuaries Coasts. 43, 1377–1393. doi: 10.1007/s12237-020-00720-5

CrossRef Full Text | Google Scholar

Hemminga M., Duarte C. (2000). Seagrass Ecology (Cambridge: Cambridge University Press).

Google Scholar

Hemminga M., Harrison P. G., Lent F. V. (1991). The balance of nutrient losses and gains in seagrass meadows. Mar. Ecol. Prog. Series. 71, 85–96. doi: 10.3354/meps071085

CrossRef Full Text | Google Scholar

Huang J. P., Fang C., Ma X., Wang L., Yang J., Luo J., et al. (2019). Tropane alkaloids biosynthesis involves an unusual type III polyketide synthase and non-enzymatic condensation. Nat. Commun. 10, 4036. doi: 10.1038/s41467-019-11987-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Invers O., Pérez M., Romero J. (2002). Seasonal nitrogen speciation in temperate seagrass Posidonia oceanica (L.) Delile. J. Exp. Mar. Biol. Ecol. 273, 219–240. doi: 10.1016/S0022-0981(02)00167-3

CrossRef Full Text | Google Scholar

IPCC. (2013). “Summary for policymakers,” in Climate Change. 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Eds. Stocker T. F., Qin D., Plattner G.-K., Tignor M., Allen S. K., Boschung J., Nauels A., Xia Y., Bex V., Midgley P. M. (Cambridge, UK and New York, NY: Cambridge University Press), 1–32.

Google Scholar

Janhs P., Latowski D., Strzalka K. (2009). Mechanism and regulation of the violaxanthin cycle: The role of antenna proteins and membrane lipids. Biochim. Biophys. Acta 1787, 3–14. doi: 10.1016/j.bbabio.2008.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang Z. J., Huang X. P., Zhang J. P. (2010). Effects of CO₂ enrichment on photosynthesis, growth, and biochemical composition of seagrass Thalassia hemprichii (Ehrenb.) Aschers. J. Integr. Plant Biol. 52 (10), 904–913. doi: 10.1111/j.1744-7909.2010.00991.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanehisa M., Goto S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28 (1), 27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch E. W., Ackerman J., Verduin J., Keulen M. v. (2007). “Fluid dynamics in seagrass ecology-from molecules to ecosystems,” in Seagrasses: Biology, Ecology and Conservation. (Dordecht: Springer), 193–225. doi: 10.1007/978-1-4020-2983-7_8

CrossRef Full Text | Google Scholar

Kolde R. (2015). Package ‘pheatmap’. R Package 1 (7), 790.

Google Scholar

Lee S. Y., Oh J. H., Choi C. I., Suh Y., Mukai H. (2005). Leaf growth and population dynamics of intertidal Zostera japonica on the western coast of Korea. Aquat. Botany. 83, 263–280. doi: 10.1016/j.aquabot.2005.07.004

CrossRef Full Text | Google Scholar

Lemos R. T., Pires H. E. (2004). The upwelling regime off the west Portuguese coast 1941-2000. Int. J. Climatology. 24, 511–524. doi: 10.1002/joc.1009

CrossRef Full Text | Google Scholar

Li B., Dewey C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323. doi: 10.1186/1471-2105-12-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin Y. X., Lin Y. F., Chen Y. H., Wang H., Shi J., Lin H. T. (2016). Hydrogen peroxide induced changes in energy status and respiration metabolism of harvested longan fruit in relation to pericarp browning. J. Agric. Food Chem. 64 (22), 4627–4632. doi: 10.1021/acs.jafc.6b01430

PubMed Abstract | CrossRef Full Text | Google Scholar

Love M. I., Huber W., Anders S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 (12), 550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Marbà N., Cebrián J., Enríquez S., Duarte C. M. (1996). Growth patterns of western Mediterranean seagrasses: speciesspecific responses to seasonal forcing. Mar. Ecol. Prog. Series. 133 (1-3), 203–215. doi: 10.3354/meps133203

CrossRef Full Text | Google Scholar

Martínez-Crego B., Olivé I., Santos R. (2014). CO₂ and nutrient driven changes across multiple levels of organization in Zostera noltii ecosystems. Biogeosciences 11, 7237–7249. doi: 10.5194/bg-11-7237-2014

CrossRef Full Text | Google Scholar

Mcleod E., Chmura G. L., Bouillon S., Salm R., Björk M., Duarte C. M., et al. (2011). A blueprint for blue carbon: Toward an improved understanding of the role of vegetated coastal habitats in sequestering CO₂. Front. Ecol. Environ. 9 (10), 552–560. doi: 10.1890/110004

CrossRef Full Text | Google Scholar

Nagelkerken I., Munday P. (2015). Animal behaviour shapes the ecological effects of ocean acidification and warming: moving from individual to community-level responses. Glob Chang Biol. 22, 974–989. doi: 10.1111/gcb.13167

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen H. M., Ralph P. J., Marín-Guirao L., Pernice M., Procaccini G. (2021). Seagrasses in an era of ocean warming: a review. Biol. Rev. Camb Philos. Soc 96 (5), 2009–2030. doi: 10.1111/brv.12736

PubMed Abstract | CrossRef Full Text | Google Scholar

Olsen J., Rouzé P., Verhelst B., Lin Y. C., Bayer T., Collen J., et al. (2016). The genome of the seagrass Zostera marina reveals angiosperm adaptation to the sea. Nature 530, 331–335. doi: 10.1038/nature16548

PubMed Abstract | CrossRef Full Text | Google Scholar

Orth R. J., Carruthers T. J. B., Dennison W. C., Duarte C. M., Fourqurean J. W., Heck K. L., et al. (2006a). A global crisis for seagrass ecosystems. BioScience 56 (12), 987–996. doi: 10.1641/0006-3568(2006)56[987:AGCFSE]2.0.CO;2

CrossRef Full Text | Google Scholar

Orth R. J., Harwell M. C., Inglis G. J. (2006b). “Chapter 5 ecology of seagrass seeds and seagrass dispersal processes,” in Seagrasses: Biology, Ecology and Conservation (Dordrecht: Springer), 111–133.

Google Scholar

Orth R. J. W., Moore K. A. (1986). Seasonal and year-to-year variation in the growth of Zostera marina L. (eelgrass) in the lower Chesapeake Bay. Aquat. Botany. 24, 335–341. doi: 10.1016/0304-3770(86)90100-2

CrossRef Full Text | Google Scholar

Parmesan C., Yohe G. (2003). A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42. doi: 10.1038/nature01286

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereda-Briones L., Tomas F., Terrados J. (2018). Field transplantation of seagrass (Posidonia oceanica) seedlings: Effects of invasive algae and nutrients. Mar. pollut. Bull. 134, 160–165. doi: 10.1016/j.marpolbul.2017.09.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Rao G., Sui J., Zhang J. (2016). Metabolomics reveals significant variations in metabolites and correlations regarding the maturation of walnuts (Juglans regia L. ). Biol. Open 5 (6), 829–836. doi: 10.1242/bio.017863

PubMed Abstract | CrossRef Full Text | Google Scholar

Richard W. R., Viva F. B., NOAA CDR Program. (2008). NOAA Optimum Interpolation 1/4 Degree Daily Sea Surface Temperature (OISST) Analysis, Version 2. [indicate subset used]. (NOAA National Centers for Environmental Information). doi: 10.7289/V5SQ8XB5

CrossRef Full Text | Google Scholar

Rischer H., Orešič M., Seppänen-Laakso T., Katajamaa M., Lammertyn F., Ardiles-Diaz W., et al. (2006). Gene-to-metabolite networks for terpenoid indole alkaloid biosynthesis in Catharanthus roseus cells. Proc. Natl. Acad. Sci. 103, 5614–5619. doi: 10.1073/pnas.0601027103

CrossRef Full Text | Google Scholar

Robinson M. D., McCarthy D. J., Smyth G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26 (1), 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Russell B. D., Connell S. D., Uthicke S., Muehllehner N., Fabricius K. E., Hall-Spencer J. M. (2013). Future seagrass beds: Can increased productivity lead to increased carbon storage? Mar. pollut. Bulletin. 73, 463–469. doi: 10.1016/j.marpolbul.2013.01.031

CrossRef Full Text | Google Scholar

Saccenti E., Hoefsloot H. C. J., Smilde A. K., Westerhuis J. A., Hendriks M. W. B. (2014). Reflections on univariate and multivariate analysis of metabolomics data. Metabolomics 10, 361–374. doi: 10.1007/s11306-013-0598-6

CrossRef Full Text | Google Scholar

Sand-Jensen K. (1975). Biomass, net production and growth dynamics in an eelgrass (Zostera marina L.) population in vellerup vig, Denmark. Ophelia 14 (1-2), 185–201. doi: 10.1080/00785236.1975.10422501

CrossRef Full Text | Google Scholar

Seddon S. (2004). Going with the flow: facilitating seagrass rehabilitation. Ecol. Manage. Restor. 5, 167–176. doi: 10.1111/j.1442-8903.2004.00205.x

CrossRef Full Text | Google Scholar

Silva J., Barrote I., Costa M. M., Albano S., Santos R. (2013). Physiological responses of Zostera marina and Cymodocea nodosa to light-limitation stress. PLoS One 8 (11), e81058. doi: 10.1371/journal.pone.0081058

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuya F., Cacabelos E., Duarte P., Jacinto D., Castro J. J., Silva T., et al. (2012). Patterns of landscape and assemblage structure along a latitudinal gradient in ocean climate. Mar. Ecol. Prog. Series. 466, 9–19. doi: 10.3354/meps09941

CrossRef Full Text | Google Scholar

Van Colen C., Ong E. Z., Briffa M., Wethey D. S., Abatih E., Moens T., et al. (2020). Clam feeding plasticity reduces herbivore vulnerability to ocean warming and acidification. Nat. Clim Change 10, 162–166. doi: 10.1038/s41558-019-0679-2

CrossRef Full Text | Google Scholar

Wang P., Laws E., Wang Y., Chen J., Song X., Huang R., et al. (2022). Elevated pCO₂ changes community structure and function by affecting phytoplankton group-specific mortality. Mar. pollut. Bull. 175, 113362. doi: 10.1016/j.marpolbul.2022.113362

PubMed Abstract | CrossRef Full Text | Google Scholar

Want E. J., Masson P., Michopoulos F., Wilson I. D., Theodoridis G., Plumb R. S., et al. (2013). Global metabolic profiling of animal and human tissues via UPLC-MS. Nat. Protoc. 8 (1), 17–32. doi: 10.1038/nprot.2012.135

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe M., Nakaoka M., Mukai H. (2005). Seasonal variation in vegetative growth and production of the endemic Japanese seagrass Zostera asiatica: A comparison with sympatric Zostera marina. Botanica Marina. 48 (4), 266–273. doi: 10.1515/BOT.2005.036

CrossRef Full Text | Google Scholar

Watson S. A., Allan B. J. M., Mcqueen D. E., Nicol S., Parsons D. M., Pether S. M. J., et al. (2018). Ocean warming has a greater effect than acidification on the early life history development and swimming performance of a large circumglobal pelagic fish. Glob Chang Biol. 24, 4368–4385. doi: 10.1111/gcb.14290

PubMed Abstract | CrossRef Full Text | Google Scholar

Waycott M., Duarte C. M., Carruthers T. J. B., Orth R. J. (2009). Accelerating loss of seagrasses across the globe threatens coastal ecosystems. Proc. Natl. Acad. Sci. U.S.A. 106, 12377–12381. doi: 10.1073/pnas.0905620106

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei T., Simko V., Levy M., Xie Y., Jin Y., Zemla J. (2017). Package ‘corrplot’. Statistician 56, 316–324.

Google Scholar

Xu M., Zhou Y., Song X., Zhang Y., Zhang H. (2019). The distribution of large floating seagrass (Zostera marina) aggregations in northern temperate zones of Bohai Bay in the Bohai Sea, China. PLoS One 14 (3), e0201574. doi: 10.1371/journal.pone.0201574

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamada Y., Sato F. (2013). Transcription factors in alkaloid biosynthesis. Int. Rev. Cell Mol. Biol. 305, 339–382. doi: 10.1016/B978-0-12-407695-2.00008-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan W., Wang Z., Pei Y., Zhou B. (2023). How does ocean acidification affect Zostera marina during a marine heatwave? Mar. Pollut. Bull. 194 (Pt B), 115394.

PubMed Abstract | Google Scholar

Yoon D., Kim Y. J., Lee W. K., Choi B. R., Oh S. M., Lee Y. S., et al. (2020). Metabolic changes in serum metabolome of beagle dogs fed black ginseng. Metabolites 10 (12), 517. doi: 10.3390/metabo10120517

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang S., Lin H., Lin Y., Lin Y., Hung Y. C., Chen Y., et al. (2017). Energy status regulates disease development and respiratory metabolism of Lasiodiplodia theobromae (Pat.) Griff. & Maubl.-infected longan fruit. Food Chem. 231, 238–246. doi: 10.1016/j.foodchem.2017.03.132

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong J., Guo Y., Liang Z., Huang Q., Lu H., Pan J., et al. (2021). Adaptation of a marine diatom to ocean acidification and warming reveals constraints and trade-offs. Sci. Total Environ. 771, 145167. doi: 10.1016/j.scitotenv.2021.145167

PubMed Abstract | CrossRef Full Text | Google Scholar