Comparative gene co-expression networks show enrichment of brassinosteroid and vitamin B processes in a seagrass under simulated ocean warming and extreme climatic events

Abstract

Introduction:

Ocean warming combined with extreme climatic events, such as marine heatwaves and flash flooding events, threaten seagrasses globally. How seagrasses cope with these challenges is uncertain, particularly for range-edge populations of species such as Posidonia australis in Shark Bay, Western Australia. Analyzing gene expression while manipulating multiple stressors provides insight into the genetic response and resilience of seagrasses to climate change. We conducted a gene expression study on a polyploid clone of P. australis during an 18-week mesocosm experiment to assess the responses to single and combined future climate change-associated stressors.

Methods:

Plants were exposed to (1) future ocean warming temperature (baseline +1.5°C) followed by a simulated marine heat wave (baseline +5.5°C), (2) light deprivation simulating observed marine heatwave driven turbidity (95% shade) at baseline temperatures, or (3) both stressors simultaneously. Basal leaf meristems were sampled for gene expression analysis using RNA-seq at four time points during the experiment. Weighted gene co-expression network analysis, GO term enrichment, and KEGG pathway enrichment analyses were used to identify stress responses.

Results:

Shaded plants showed specific gene enrichment for shade avoidance (programmed cell death) after three weeks of stress, and before any heated tanks showed a specific heat response. Shaded plants were positively correlated with programmed cell death and stress-related processes at the end of the experiment. Once ocean warming temperatures (+1.5°C) were in effect, gene enrichment for heat stress (e.g., ROS scavenging and polyamine metabolism) was present. Vitamin B processes, RNA polymerase II processes. and light-related meristematic phase changes were expressed with the addition of simulated MHW. Heated plants showed meristematic growth signatures as well as trehalose and salicylic acid metabolism. Brassinosteroid-related processes were significantly enriched in all stressor treatments at all time points, except for the isolated heat-stressed plants three weeks after stressor initiation.

Discussion:

Gene expression responses to the interaction between heat waves and turbidity-induced light reduction support the observed geographical scale mortality in seagrasses observed for P. australis in Shark Bay, suggesting that even this giant polyploid clone will be negatively impacted by more extreme climate change projections.

Introduction

Ecosystems are affected by global climate change (e.g., Scheffers et al., 2016; Hughes et al., 2017). Some of these changes are clearly visible and have received much public attention, such as coral bleaching events (Hughes et al., 2017), whereas other changes are more insidious and hidden. Global temperatures are projected to continue to rise throughout the 21st century with extreme climatic events becoming more common. Marine heat waves (MHWs) have been predicted to increase in frequency and intensity (Hobday et al., 2018; Oliver et al., 2018). Such extreme events can cause widespread mortality and collapse of ecosystems and their services (Wernberg et al., 2016; Pecl et al., 2017). Identifying how individuals respond to climate-induced stressors that negatively affect the resilience of species is fundamental to managing the impacts of climate change (Hoffmann and Sgrò, 2011; Lankau et al., 2011).

Environmental stressors affect plant physiology and fitness through changes in morphology, reproduction, survival, growth rates, metabolic rates, and gene expression (Hoffmann and Sgrò, 2011; Zandalinas et al., 2022). Rising ocean temperatures and reduced light availability, for example, caused by flash flooding and subsequent sediment resuspension in the water column, are two key stressors for marine plants (Duarte et al., 2018). Each stressor has an effect in isolation; however, in combination, they can have additive or synergistic effects (Brook et al., 2008; Jung et al., 2023). Interactions among stressors can push seagrass populations above their functional threshold, where range contractions and extinctions are likely to occur (York et al., 2013; Lefcheck et al., 2017; Kendrick et al., 2019). Gene expression studies have provided useful information on climate stressors in isolation, but only recently have stressors been tested in combination (Zhang and Sonnewald, 2017; Tsioli et al., 2022; Santillán-Sarmiento et al., 2023).

Genomic tools can be combined with ex situ mesocosm experiments to increase our understanding of the ability of a species to respond to climate change (Stillman and Armstrong, 2015). Transcriptomic studies have indicated that gene regulation plays a considerable role in phenotypic plasticity and acclimation (King et al., 2018; Rivera et al., 2021). Changes in the expression levels of candidate genes have been reported in mesocosm experiments using seagrasses, such as those associated with single stressors (Franssen et al., 2011; Franssen et al., 2014) and interactions among multiple stressors (Marín-Guirao et al., 2016; Malandrakis et al., 2017; Zhang et al., 2023). However, with the introduction of RNA-seq capabilities, measuring global transcriptomic changes provides an opportunity to gain a broader understanding of the gene networks that influence stress responses.

A polyploid Posidonia australis clone was recently discovered in the metahaline waters of Shark Bay (Edgeloe et al., 2022), persists across 200 km², and has a wide environmental gradient (15°C–27°C, 35 PSU to 49 PSU, high light exposure, 1,500 µmol m s⁻¹ to 3,000 µmol m s⁻¹; Edgeloe et al., 2022). Climate change has been identified as the single largest threat to the persistence of temperate seagrasses in Shark Bay (Kendrick et al., 2019), with ~36% of meadows affected during the last extreme marine heatwave (MHW) in the Austral summer of 2010–2011 (Strydom et al., 2020).

We conducted an 18-week mesocosm experiment using the polyploid P. australis clone from Shark Bay to assess temporal gene expression responses to isolated and combined effects of elevated heat and light deprivation stressors in the context of an MHW under future climate change projections. We addressed the following research questions: (1) What are the effects of projected warming with climate change in combination with temperature increases associated with a MHW? (2) What are the effects of extreme reductions in light associated with projected climate change and MHW increases in storm events and flooding-induced turbidity? (3) What are the additive and/or synergistic effects of both projected increases in temperature and decreased light associated with climate-induced ocean warming and the increased frequency of marine heatwaves? We used RNA-seq and weighted gene co-expression network analysis to identify stressor-responsive genes that were significantly positively correlated with low light (shade), elevated heat (heat), and a combination of the two (heat and shade) throughout the experiment.

Methods

Plant collections

Ramets of the polyploid P. australis clone were harvested from Middle Bluff (25.823980°S, 113.464010°E) in the western gulf of Gathaagudu (in the local Malgana language), Shark Bay, northwest Western Australia, in March 2020. Approximately 350 ramets were collected via SCUBA from growing edges of the meadow at an average depth of 1.5 m. The in situ temperature was 26.1°C ± 0.3°C (± standard error) and salinity was 40.2 PSU ± 0.2 PSU (Practical Salinity Units) at the time of collection. The details of transportation and planting are described in Jung et al. (2023).

Mesocosm experiment

A mesocosm experiment was conducted over 18 weeks (Figure 1). A detailed description of the mesocosm setup can be found in Jung et al. (2023). Briefly, 12 L × 1,800 L self-contained, recirculating, reinforced tanks were in a temperature-controlled glasshouse with ambient light and natural seawater. Twenty-seven pots of single P. australis ramets containing one apical shoot and up to four mature shoots were placed in each tank for a 10-week acclimation period to ensure sufficient time for plant adjustment to the mesocosm. The water temperature and salinity were maintained at long-term averages for the Middle Bluff of 26°C (Wijffels et al., 2018) and 42 PSU ± 3 PSU.

Figure 1

After the 10-week acclimation period, three tanks each were randomly assigned to either (A) ‘Control’ at a baseline temperature of 26°C and ambient light levels or one of three treatments: (B) a high temperature stressor designated ‘Heat’, (C) a low light stressor designated ‘Shade’, or (D) combined interaction between high temperature and low light stressors designated ‘Heat + Shade’ (Figure 1). Water temperature in the ‘Heat’ and ‘Heat + Shade’ treatments was elevated +1.5°C over two days to reach a future warming temperature simulating the projected minimum end-of-century temperature rise, equivalent to 27.5°C in Shark Bay (Geraldi et al., 2020; IPCC, 2021). The six ‘Shade’ and ‘Heat + Shade’ treatment tanks were overlaid with a double layer of high-density shade cloth to achieve 95% light deprivation simulating high turbidity levels observed in 2011 in Shark Bay, (Fraser et al., 2014). The shade cloth was not removed from ‘Shade’ or ‘Heat + Shade’ treatment tanks until the end of the experiment. Following the 3-week future warming period (baseline +1.5°C), both the ‘Heat’ and ‘Heat + Shade’ tanks were then increased by +4.0°C in +0.8°C increments over five days, simulating the +4°C increase in water temperatures experienced in Shark Bay during the 2011 heatwave (Fraser et al., 2014; Kendrick et al., 2019). After two weeks of MHW, temperatures were decreased in ‘Heat’ and ‘Heat + Shade’ tanks by −0.8°C increments over another five days, back to future warming temperatures (+1.5°C above baseline control) and held for one week to simulate a recovery period. Temperature, light, and salinity levels were monitored at the seagrass leaf canopy height throughout the experiment, as reported in the Supporting Information in Jung et al. (2023).

Sampling for gene expression response

Three pots containing a single P. australis ramet with one apical shoot and up to four mature shoots were randomly harvested per tank at four time points to assess the gene expression responses. The four sampling efforts were as follows: T0: after ten weeks of acclimation at a baseline 26°C, T1: after three weeks of future warming temperature (baseline +1.5°C) and/or three weeks of shade, T2: after one week of simulated MHW (baseline +5.5°C) and/or five weeks of shade, and T3: after one week of MHW recovery under future warming temperature (baseline +1.5°C) and/or eight weeks of shade (Figure 1). Approximately 2 cm of basal leaf meristem tissue was collected from a single shoot, corresponding to the leaves used in Jung et al. (2023) (youngest fully developed leaf from the apical shoot), and immediately frozen in liquid nitrogen prior to storage at −80°C.

RNA extraction, sequencing, and assembly

Approximately 100 mg of basal leaf meristem tissue was weighed for each sample (three samples from each treatment from four time points, n = 48 samples) and ground in a 2 mL screw-cap tube with six 3 mm spherical YSZ grinding media beads using a 2010 Geno/Grinder^® (SPEX Sample Prep, Metuchen, NJ, USA). Samples were placed in prechilled (-80°C) metal tube blocks while grinding and left for 1 min at 1,300 rpm or pulsed for 30 s until sufficiently ground into powder. The tube block was refreshed with liquid nitrogen to ensure that samples remained chilled. A Spectrum™ Plant Total RNA kit (Sigma-Aldrich^® St Louis, MO, USA) was used to extract total RNA from the samples, which was then quantified and treated according to the protocols described by Booth et al. (2022). RNA was sent for poly(A)-selected library preparation for mRNA sequencing via an Illumina^® HiSeq^® sequencing platform (short reads of 2 bp × 150 bp at 30× coverage) by GENEWIZ^® (Suzhou, China).

Raw fastQ sequence reads were initially assessed for quality using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and were mapped to the P. australis genome assembly (GCA_024699735.1; Bayer et al., 2022 preprint). Mapping was performed using HISAT2 v2.0.5 (http://daehwankimlab.github.io/hisat2/) (Kim et al., 2019) and samtools v1.9 (http://samtools.sourceforge.net/) (Li et al., 2009). The mapped reads were assembled into transcripts using the StringTie 2.1.1 (https://ccb.jhu.edu/software/stringtie/; Pertea et al., 2015). Gffcompare (https://ccb.jhu.edu/software/stringtie/gffcompare.shtml; Pertea and Pertea, 2020) was used to compare StringTie transcripts with known transcripts.

Temporal gene co-expression networks

Weighted Gene Co-expression Network Analysis (WGCNA) (Langfelder and Horvath, 2008) was carried out using the R packages WGCNA v1.71 and DESeq2 v1.34.0 (Love et al., 2014). WGCNA creates networks of gene expression patterns based on count data generated using the publicly available Python script prepDE.py3 (http://ccb.jhu.edu/software/stringtie/dl/prepDE.py3) and custom phenotypic metadata using Microsoft Excel. Gene counts were filtered to include only genes with counts ≥15 in at least 75% of the samples, then normalized using the “vst” (variance stabilizing transformation) function of the DESeq2 package. Genes were assigned to modules based on similar expression patterns that were correlated with the experimental treatments. Four co-expression networks were created, one for each time point (T0, T1, T2, and T3) where each experimental treatment was described: 'Heat', 'Shade', 'Heat + Shade', and ‘Control’. Principal Component Analysis (PCA) was used to visually inspect the data via the R package ggplot2 v3.3.6 (Wickham and Wickham, 2016). All networks were “signed” networks where the soft-thresholding power (the power that correlations among genes are raised to reduce noise) was set to a conservative value of 18 for all networks, as recommended per WGCNA guidelines, since no soft-threshold values were appropriate (decided by simultaneously having the lowest value to meet our 0.8 threshold for scale free topology and sufficiently close to zero mean connectivity score). These settings produced conservative networks that lack a scale-free topology fit owing to interesting biological variables (groups of samples). An exception was made for the T2 network, which fit the topology and connectivity models at a value of 20 (Supplementary Figures 1A–D). Once the power was chosen we then used the “BlockwiseConsensusModules” function of WGCNA to detect modules and construct the networks with default values except: “TOMtype” chosen was “signed,” “maxBlockSize was 5,000 and the soft threshold given as above. All heatmaps were generated using the “CorLevelPlot” CorLevelPlot R package (Blighe, 2018) which specifies the Pearson correlation method used to calculate the correlation coefficients between module eigengenes and experimental treatments. The blue to red scale was used to designate negative to positive correlations.

Gene ontology enrichment and KEGG pathway enrichment

Genes were extracted from modules that were shown to have a significant correlation with each given treatment in the module-treatment heatmaps for all networks. These genes were then annotated for Gene Ontology (GO) terms using PANNZER2 (http://ekhidna2.biocenter.helsinki.fi/sanspanz/; Koskinen et al., 2015) and KEGG ID via KAAS (https://www.genome.jp/tools/kaas/; Moriya et al., 2007). The R package TopGO v2.46.0 (Alexa and Rahnenfuhrer, 2021) was used to conduct GO enrichment with a node size of five for genes in their corresponding modules, to produce a list of gene ontologies (classifications of gene functional annotations) found significantly enriched in a module/treatment of interest relative to all samples. The enriched GO terms and their P-values were input to REVIGO (http://revigo.irb.hr/; Supek et al., 2011) with a SimRel cut-off value of 0.7 against the Uniprot database (https://www.uniprot.org) to generate a condensed list of enriched GO terms. The R package clusterProfiler v.4.2.2 (Wu et al., 2021) was then used to carry out pathway enrichment on the same modules of genes using Benjamini–Hochberg P-value adjustment. This process, similar to GO enrichment, produces a list of molecular pathways (classified by KEGG database pathway annotation) found to be significantly enriched in a module/treatment of interest relative to all samples.

Comparative presence of metabolite activity

The presence of genes involved in the biological processes of metabolites of interest from Jung et al. (2023) was searched in the transcriptome generated in this study, in accordance with the sister experiment conducted by Jung et al. (2023). A custom program was made using R to search the terms: “glucose,” “trehalose,” “sucrose,” “glucopyranose,” “polyol,” “ribitol,” “dulcitol,” “erythritol,” “ribonic,” “hexaconic,” “L-pyrolglutamate,” “L-aspartate,” “l-threonine,” “serine,” “myo-inositol,” “shikimic,” “shikimate,” “isoleucine,” “alanine,” “valine,” “phenylalanine,” “putrescine,” “rhamnose,” “glycerol,” “proline,” “threonic,” “gluconic,” “aspartate,” “citric,” “fructose,” “amino acid,” “sugar,” and “methionine” to capture all relevant metabolite-related biological processes in the enriched GO terms of modules significantly positively correlated to any treatment across three exposure time points (T1–T3). The matched terms were then filtered for relevance and accuracy.

Results

Transcriptome sequencing and assembly

Extracted RNA was sequenced at a depth of 20 million reads, resulting in a total of 2,069,079,918 raw 150 bp paired-end reads across all samples (Supplementary Table 1). FastQC analysis showed that all samples contained mean Phred scores >36, indicating that these were high-quality reads with no need for trimming and filtering owing to redundancy in downstream programs. On average, 74% of the reads for each sample were mapped to the P. australis genome (Bayer et al., 2022 preprint; Supplementary Table 1), resulting in 42,568 non-redundant genes. The assembly statistics of the P. australis genome used in this study contained a total size of 1,215 Mbp, an N50 score of 9,415, and 258,843 contigs longer than 1 kbp (Bayer et al., 2022 preprint). The PCA plot showed no clear clustering among the samples following a 10-week acclimation period at 26°C (T0; Figure 2A). There was also no significant positive correlation between any module and any treatment by WGCNA (T0, Figure 3A). In addition, there were no significantly enriched pathways at this time point. Approximately 28% (11,972 of 42,568) of total genes detected were significantly positively correlated to the three stressor treatments (‘Heat’, ‘Shade’, and ‘Heat + Shade’) across the sampled time points (T1–T3). Our interest was on stressor-related biological processes and the translational evidence of specific metabolites found in relation to ‘Heat’, ‘Shade’, and ‘Heat + Shade’ treatments. Therefore, only the significantly positive correlations between modules and treatments were further explored.

Figure 2

Figure 3

Temporal change within single stressor treatments

The overall effect of the ‘Heat’ treatment was not as apparent early in the experiment (baseline +1.5°C, T1) as the ‘Shade’ treatment, as shown in comparison of T0 and T1 PCA plots (Figures 2A, B) where ‘Heat’ and ‘Control’ tanks remain loosely clustered in T0 (PC1 78% variance, PC2 6% variance) and T1 (PC1 78% variance, PC2 6% variance), whereas both ‘Shade’ and ‘Heat + Shade’ tanks cluster more tightly than the ‘Control’ and ‘Heat’ tanks along PC1 and PC2 of both plots. At T2 during the peak of the MHW (baseline +5.5°C), both ‘Heat’ and ‘Heat + Shade’ tanks clustered together of PC1 (69% variance) and independent of each other along PC2 (12% variance), though one replicate each of ‘Shade’ joined one replicate of ‘Control’ remaining together from T1 to T2 (Figures 2B, C). By the recovery period at T3, both the ‘Heat + Shade’ and ‘Heat’ replicates clustered the most tightly together along both PC1 (78% variance) and PC2 (12% variance) with two replicates of the ‘Control’ tanks and one from ‘Shade.’ A single ‘Control’ replicate remained apart from all other samples along PC1, however, now two replicates of the ‘Shade’ tanks also appear apart from all other samples along PC1, and independent of one another along PC2.

Genes modules were unique to a treatment or shared between treatments within a given time point (Figures 3A–D; Supplementary Table 2). Unique gene modules were observed for each treatment during the recovery period (T3), with some shared among all treatments (Figure 3D; Supplementary Table 2). We considered the repeated presence of significantly positively correlated with biological processes over multiple time points as a sustained effect of the stressor during the experiment. Modules found significantly positively correlated to one of the individual treatments and ‘Heat + Shade’ within the same time point were reported in the ‘Heat + Shade and Shared Modules’ section below. No modules were found significantly positively correlated in both individual treatments (‘Heat’ or ‘Shade’) and not in the ‘Heat + Shade’ treatment. Genes were considered a generalized stress response if the module was significantly positively correlated with all treatments.

‘Heat’ modules

No modules in the WGCNA module-treatment heatmap exhibited significantly positive correlations with ‘Heat’ after three weeks at future warming temperature (baseline +1.5°C, T1; Figure 3B). After three weeks of future warming plus one week of simulated MHW, three modules were observed to be significantly positively correlated with ‘Heat’ (baseline + future warming +4°C, T2): “gray”: 0.91, P <0.05; “magenta”: 0.88, P <0.05 and “darkgray”: 0.96, P < 0.01 (Figure 3C). These three modules were also positively correlated with the ‘Heat + Shade’ treatment and are reported in the ‘Heat + Shade and Modules Shared Within Time point’ section. After one week of MHW recovery, where temperatures were reduced to baseline +1.5°C at T3, one module was uniquely significantly positively correlated in the ‘Heat’ treatment (“sienna3”: 0.84, P <0.05; Figure 3D). This is the only module unique to the ‘Heat’ treatment and contains genes significantly enriched in iron homeostasis and iron transport, as well as being significantly enriched in the plant hormone signal transduction KEGG pathway (Supplementary Table 2).

‘Shade’ modules

At three weeks of exposure to 95% shade stress (T1), a single module was significantly positively correlated to the ‘Shade’ treatment (“darkorange”: 0.87, P <0.05; Figure 3B). The genes in this module showed enriched biological processes involved in brassinosteroid signaling and growth regulation (Supplementary Table 2). These terms are largely driven by a single, novel gene. KEGG pathway enrichment analysis identified the aminoacyl-tRNA biosynthesis pathway as significantly enriched within this module (Supplementary Table 2). At T2, after five weeks of shading, a single module was uniquely significantly positively correlated (“royalblue”: 0.95, P <0.01; Figure 3C). The genes in this module primarily showed enrichment for leaf senescence and related programmed cell death (Supplementary Table 2). By T3 (after eight weeks of shade), the single module unique to the ‘Shade’ treatment (“greenyellow”: 0.84, P <0.05; Figure 3D) was significantly enriched in terms involved in accumulation of chloroplasts, as well as mitochondrial respiratory chain in relation to shaded conditions.

‘Heat + shade’ and modules shared within time point

At T1, three weeks into the experiment at +1.5°C and 95% shading, two modules were significantly positively correlated to both the ‘Shade’ (“brown”: 0.95, P <0.01; “steelblue”: 0.92, P <0.05) and ‘Heat + Shade’ treatments (“brown”: 0.92, P <0.05; “steelblue”: 0.85, P <0.05; Figure 3B). No specific heat-related processes were significantly enriched in either module. The “brown” module was significantly enriched in a broad range of GO terms involved in programmed cell death, TCA cycle, and shade avoidance, along with the spliceosome KEGG pathway also found enriched in this module (Supplementary Table 2). Additionally, the “steelblue” module did not contain genes significantly enriched in an obvious over all biological process or pathway, with prominent terms such as “telomere maintenance,” “DNA damage response,” and “mitochondrial transmembrane transport”.

At T2 (one week of MHW and four weeks of shade exposure), three modules were significantly positively correlated to the ‘Heat + Shade’ treatment’ (“pink”: 0.93, P < 0.01; “turquoise”:0.83, P <0.05; “yellow”:0.84, P <0.05; Figure 3C). The “pink” module is significantly enriched in biological processes involved in RNA polymerase II transcription and translation, as well as light-related processes and meristem vegetative to reproductive phase transition with terms such as “cellular response to blue light,” “long-day photoperiodism, flowering,” “regulation of meristem growth,” “positive regulation of reproductive process,” and “vegetative to reproductive phase transition in meristem.” The “turquoise” module was found to be significantly enriched in a broad number of processes including those involved in tetrahydrofolate (vitamin B9) interconversion and purine nucleoside/nucleotide biosynthesis, as well as meristematic development and maintenance and ROS scavenging (Supplementary Table 2). The “yellow” module instead showed a broad range of stress related biological processes such as “regulation of cellular response to stress,” “response to endoplasmic reticulum stress,” “response to temperature synthesis” as well as calcium, gibberellic acid, and phosphate processes.

Two modules were found significantly positively correlated between both the ‘Heat’ (“gray”: 0.91, P <0.05; “magenta”: 0.88, P <0.05) and ‘Heat + Shade’ treatments (“gray”: 0.93, P <0.01; “magenta”: 0.92, P <0.01; Figure 3C; Supplementary Table 2). However, both the “gray” and “magenta” modules did not show a clear biological theme upon GO term enrichment. Lastly, the “darkgray” module was shown to be significantly positively correlated to all treatments (‘Heat’: 0.96, P <0.01; ‘Shade’: 0.92, P <0.01; and ‘Heat + Shade’: 0.96, P <0.01), and as such was characterized as a module representing a generalized stress response (Figure 3C). This module was found significantly enriched in the terms: “DNA topological change,” “protein glycosylation,” “brassinosteroid mediated signaling pathway,” and “response to chitin”.

At the end of a one-week recovery period at future warming temperature +1.5°C (T3) and 8 weeks of shade, two modules were significantly positively correlated to the ‘Heat + Shade’ treatment (“salmon”: 0.85, P <0.05; “darkgray”: 0.94, P <0.01; Figure 3D; Supplementary Table 2). The “salmon” module showed processes enriched in RNA transcription and translation, while other terms were involved in meristematic growth. The “darkgray” module was largely enriched in GO terms involved in the production/transport of complex sugars including carbohydrate and starch. Programmed cell death-related processes were also enriched in this module.

The “darkred” (‘Heat’: 0.88, P <0.05; ‘Shade’: 0.88, P <0.05; ‘Heat + Shade’: 0.9, P <0.05) and “magenta” (‘Heat’: 0.89, P <0.05; ‘Shade’: 0.84, P <0.05; ‘Heat + Shade’: 0.92, P <0.05) modules were both shown to be general stress response modules, being significantly positively correlated with all treatments (Figure 3D). The “darkred” module contained terms significantly enriched in processes involving brassinosteroid signaling, nitrogen, and lateral root morphogenesis. The “magenta” module instead was found to be significantly enriched in terms involved in general stress responses in plants, such as trehalose and salicylic acid metabolism, nitrogen, as well as RNA translation regulation.

Comparative presence of metabolite activity

Only seven out of 103 matched terms remained after filtering for GO enrichment relevance with P <0.05, and five were found in modules that were significantly positively correlated across all treatments (Table 1). Sampling at T1 (+1.5°C), Glucose-6-Phosphate, the first intermediate of glucose metabolism, and S-adenosylmethionine, key in the oxidative stress response, were significantly enriched (‘Shade’, ‘Heat + Shade’). With the addition of a MHW (T2, baseline +5.5°C), Glycerol-3-Phosphate (‘Heat’, ‘Heat + Shade’), the output of the Calvin Cycle and source of carbohydrates was significantly enriched. Similar trends were observed for polyol Sorbitol (‘Heat + Shade’) which is integral for the conversion of glucose to fructose, and for L-phenylalanine (‘Heat + Shade’), a protein building block that can act as a temperature stress defense compound. After one week recovery from the MHW (baseline +1.5°C), L-phenylalanine was still significantly enriched along with S-adenosylmethionine, a precursor to ethylene production and plant senescence, and Trehalose, a protectant against abiotic stressors including heat, high salinity, and others (‘Heat’, ‘Shade’, ‘Heat + Shade’).

Discussion

This study is one of the few that attempted to identify early warning responses that occur soon after stressor initiation via RNA-seq to climate change stressors in marine plants and is strengthened by the combination of gene expression and metabolomics in the same experiment. The largest observed changes in gene expression appeared under potential additive or synergistic stress conditions that simulated dramatic decreases in light (95% reduction) and extreme increases in temperature (baseline +5.5°C). The greater impact of synergistic rather than single stressors on biological systems has previously been reported in observations of climate change in marine (Oomen and Hutchings, 2017) and terrestrial (Zandalinas and Mittler, 2022) ecosystems. Our investigation of changes in expressed genes with projected future warming temperatures (baseline +1.5°C) was minimal, with predominantly generalized processes commonly seen in stress responses. We attribute this finding to the plasticity of natural gene expression responses to environmental stimuli, which is not indicative of a true stress response. However, genes associated with a stress response to additional warming by simulating a marine heat wave (baseline +5.5°C) showed an increased metabolic load. Post-MHW recovery was minimal, and a sustained stress response was observed across all treatments in trehalose (Sadak et al., 2019) and salicylic acid (El-Esawi et al., 2017) metabolism. These metabolites were also found in the regulation studied in the same experiment (Jung et al., 2023), and have been linked to the activation of antioxidant defense against reactive oxygen species (ROS) and the cellular stress-induced damage they cause, a common effect of almost any environmental stressor on plant health (Sewelam et al., 2016). While we acknowledge that no parameter other than RNA-seq was measured to corroborate a stress response, by the end of the experiment, we observed high mortality and/or signs of leaf damage across treatments, except in the control tanks. However, even with this conserved response, these results suggest a +5.5°C temperature increase would impact seagrass performance and may lead to increased seagrass morbidity and mortality.

The interaction between shade and heat had a unique stress response signature at 5.5°C above the baseline, involving genes for vitamin B9 biological processes. B-vitamins are metabolite precursors that break down during stress, leading to the inhibition of metabolic processes (Hanson et al., 2016). The role of vitamin B, specifically tetrahydrofolate (vitamin B9) and its one-carbon reactions, has received minimal attention in relation to plant stress, including in seagrasses. Vitamin B is important due to its role as metabolic co-factor precursors (Hanson et al., 2016). Vitamin B deficiency may be caused by abiotic stressors in plants, and tetrahydrofolate pathways have been reported to modulate gene expression under multiple stressors (e.g., salinity, drought, oxidative stress; Hanson et al., 2016). Differentially regulated vitamin B expression was also found during oxidative stress in Arabidopsis (Baxter et al., 2007), which is consistent with our findings. We suggest the need for further vitamin B research in seagrass stress responses, as their role may go beyond that of metabolic precursors.

A sustained response of brassinosteroid signaling-related processes was observed at all time points in almost all experimental treatments. The basal leaf meristem tissue is a hub for growth and signal transduction, and brassinosteroid signaling plays a central role in plant growth (Singh and Savaldi-Goldstein, 2015; Booth et al., 2022). However, the full context of brassinosteroid action under environmental stress is not yet fully understood (Saini et al., 2015). Brassinosteroid signaling in plants has been shown to promote plant growth, development and abiotic stress responses, in particular, they have been shown to help alleviate stress caused by low-light conditions through growth-mediated shade avoidance systems (Keuskamp et al., 2011). Similarly, brassinosteroids have been linked to higher chloroplast levels in rice, and thereby increasing heat tolerance (Pantoja-Benavides et al., 2022). Brassinosteroid metabolism was found to be upregulated in two temperate seagrasses (Nanozostera noltii and Zostera marina) which did not exhibit a strong sympatric heat response (Franssen et al., 2014). Further research is needed to understand the relative abundance and potentially protective interactions between brassinosteroids and chloroplast abundance under shading stress.

Programmed cell death and senescence responses were quite immediate, occurring in shading and the interaction between heat and shading well before any heat response. Evidence of a sustained programmed cell death/shade avoidance response in ‘Shade’ tanks is consistent with shading stress not being abated once applied. Low-light or shaded conditions can be beneficial for plant health, where short-day or shaded plants produce less ROS than under high-light conditions, which has many flow-on effects including increased resilience to thermal, drought, and biotic pathogen stressors (Roeber et al., 2021). However, in our experiment, where plants were subjected to sustained 95% shade, we suggest that programmed cell death responses are a result of circadian stress, in which detrimental ROS gradually accumulate over a prolonged period of darkness (Nitschke et al., 2016). Prolonged-dark or circadian stress studies largely involve darkness exposure, or in our case, a 95% reduction in light (shade), in the order of hours or days in established plants. In our experiment, the first time point post-stressor initiation (T1) was measured after three weeks of shading. The direct effects of prolonged darkness phenomena have been understudied, especially in monocots (Nitschke et al., 2016; Sobieszczuk-Nowicka, 2017), and should be a source of further research for other areas where turbid waters may impact seagrass meadows, such as Shark Bay (Fraser et al., 2014; Kendrick et al., 2019), or other marine environments increasingly impacted by anthropogenic disturbances and high rainfall events associated with climate change.

Transcriptomic and metabolomic results were used to verify each other to increase the reliability of enriched genes influencing the metabolome both in our study and in another recent study of the effects of heat waves on another temperate seagrass, Z. marina (Zhang et al., 2021). Gene enrichment under ocean warming and climate change-induced shading supported the metabolite regulation observed in the same experiment (Jung et al., 2023). Our study with P. australis showed that enrichment of gene expression in response to projected ocean warming and reduced light initially resulted in enrichment of genes that are associated with increased metabolic activity and oxidative stress (Glucose-6-Phosphate, S-adenosylmethionine; Amir and Hacham, 2008; Rende et al., 2019). Gene response to the addition of a marine heat wave resulted in the enrichment of genes associated with glucose metabolism and stress responses (Glycerol-3-Phosphate, Sorbitol, L-phenylalanine; Krasensky and Jonak, 2012; Yoo et al., 2013; Lai et al., 2015). No recovery from MHW was observed within one week as observed through the enrichment of stress response genes and a precursor gene for plant senescence (L-phenylalanine, trehalose, methionine; Grennan, 2007; Amir and Hacham, 2008; Yoo et al., 2013). These gene responses were consistent with the regulation of metabolites, and Jung et al. (2023) reported an increase in Glycerol-3-Phosphate, L-phenylalanine and a decrease in trehalose as a result of shade alone or in combination with MHW stress. Similarly, observations of seagrass shoot thinning in P. australis and defoliation in another temperate seagrass (Amphibolis antarctica) (Fraser et al., 2014; Kendrick et al., 2019) from the additive effects of baseline +5.5°C warming and turbidity-driven light deprivation generated in the 2010–2011 extreme MHW in Shark Bay are consistent with transcriptomic and metabolomic results from this mesocosm experiment. In combination, these stressors can result in geographical scale mortality in seagrasses (Fraser et al., 2014), where the greatest losses occur when extreme heating persists for ≥94 days (Strydom et al., 2020).

We examined gene expression in a single polyploid P. australis clone from the northern (warm) range edge of Western Australia. Individual shoots from this clone had a low range of phenotypic variation, with 28% of the transcribed genes significantly correlated to modules that correlated with the experimental stressors. It has been suggested that the clone originated as a result of historical autopolyploidy, with more recent gene loss through rediploidization, although high heterozygosity levels (~90%) are more consistent with hybridization with other Posidonia species through allopolyploidy (Edgeloe et al., 2022). In any case, our results suggest that the genome of this P. australis clone may lack vital genes or regulation among the additional gene copies needed to survive prolonged shade or combined heat and shade. This clone was showing early signs of recovery in ‘Heat’ treatment tanks by T3; however, a longer recovery time would be required to gain a better understanding of recovery capacity. Whether this example of polyploidy in a seagrass is more tolerant to abiotic stressors representative of P. australis more widely (Sinclair et al., 2023) and/or its diploid progenitors (see Edgeloe et al., 2022), as polyploidy is often thought to confer (Van de Peer et al., 2021; Heslop-Harrison et al., 2023), requires further comparative study in this context with the neighboring diploid P. australis.

In conclusion, brassinosteroids were consistently enriched throughout the heat and/or shade stress responses. We observed that this P. australis clone in Shark Bay could withstand future ocean warming temperatures (+1.5°C) but showed a significant transcriptomic response to the additional pressure of an extreme climatic event (MHW). The effect of prolonged shading was immediate and extreme, with stressor-induced damage control being a major theme in the correlated gene expression. This damage response was further exacerbated by the interaction of ‘Heat’ and ‘Shade’, inducing a novel vitamin B gene expression response. However, meaningful plant recovery was not observed at the end of this experiment. Therefore, it is unlikely that this P. australis clone would survive extreme marine heatwave conditions (Hobday et al., 2018) under future ocean warming. Further research into the role of brassinosteroids and vitamin B in relation to plant stress responses will be important for understanding the resilience of P. australis over longer recovery phases.

References

• 1

  AlexaA.RahnenfuhrerJ. (2021). TopGO: enrichment analysis for gene ontology. 2.46. 0. R Package. Available at: https://bioconductor.org/packages/topGO.

  • Google Scholar

• 2

  AmirR.HachamY. (2008). Methionine metabolism in plants. Sulfur: A Missing Link Between Soils Crops Nutr.50, 251–279. doi: 10.2134/agronmonogr50.c16

  • CrossRef
  • Google Scholar

• 3

  BaxterC. J.RedestigH.SchauerN.RepsilberD.PatilK. R.NielsenJ.et al. (2007). The metabolic response of heterotrophic Arabidopsis cells to oxidative stress. Plant Physiol.143 (1), 312–325. doi: 10.1104/pp.106.090431

  • CrossRef
  • Google Scholar

• 4

  BayerP. E.FraserM. W.MartinB. C.PetereitJ.Severn-EllisA. A.SinclairE. A.et al. (2022). Not all pathways are the same-unique adaptations to submerged environments emerge from comparative seagrass genomics. doi: 10.1101/2022.11.22.517588. Preprint.

  • CrossRef
  • Google Scholar

• 5

  BligheK. (2018). CorLevelPlot: Visualise correlation results and test significancies of these. Available at: https://github.com/kevinblighe.

  • Google Scholar

• 6

  BoothM. W.BreedM. F.KendrickG. A.BayerP. E.Severn-EllisA. A.SinclairE. A. (2022). Tissue-specific transcriptome profiles identify functional differences key to understanding whole plant response to life in variable salinity. Biol. Open11 (8), bio059147. doi: 10.1242/bio.059147

  • CrossRef
  • Google Scholar

• 7

  BrookB. W.SodhiN. S.BradshawC. J. (2008). Synergies among extinction drivers under global change. Trends Ecol. Evol.23 (8), 453–460. doi: 10.1016/j.tree.2008.03.011

  • CrossRef
  • Google Scholar

• 8

  DuarteB.MartinsI.RosaR.MatosA. R.RoledaM. Y.ReuschT. B.et al. (2018). Climate change impacts on seagrass meadows and macroalgal forests: an integrative perspective on acclimation and adaptation potential. Front. Mar. Sci.5, 190. doi: 10.3389/fmars.2018.00190

  • CrossRef
  • Google Scholar

• 9

  EdgeloeJ. M.Severn-EllisA. A.BayerP. E.MehraviS.BreedM. F.KraussS. L.et al. (2022). Extensive polyploid clonality was a successful strategy for seagrass to expand into a newly submerged environment. Proc. R. Soc. B289 (1976), 20220538. doi: 10.1098/rspb.2022.0538

  • CrossRef
  • Google Scholar

• 10

  El-EsawiM. A.ElansaryH. O.El-ShanhoreyN. A.Abdel-HamidA. M.AliH. M.ElshikhM. S. (2017). Salicylic acid-regulated antioxidant mechanisms and gene expression enhance rosemary performance under saline conditions. Front. Physiol.8, 716. doi: 10.3389/fphys.2017.00716

  • CrossRef
  • Google Scholar

• 11

  FranssenS. U.GuJ.BergmannN.WintersG.KlostermeierU. C.RosenstielP.et al. (2011). Transcriptomic resilience to global warming in the seagrass Zostera marina, a marine foundation species. Proc. Natl. Acad. Sci.108 (48), 19276–19281. doi: 10.1073/pnas.1107680108

  • CrossRef
  • Google Scholar

• 12

  FranssenS. U.GuJ.WintersG.HuylmansA. K.WienpahlI.SparwelM.et al. (2014). Genome-wide transcriptomic responses of the seagrasses Zostera marina and Nanozostera noltii under a simulated heatwave confirm functional types. Mar. Genomics15, 65–73. doi: 10.1016/j.margen.2014.03.004

  • CrossRef
  • Google Scholar

• 13

  FraserM. W.KendrickG. A.StattonJ.HoveyR. K.Zavala-PerezA.WalkerD. I. (2014). Extreme climate events lower resilience of foundation seagrass at edge of biogeographical range. J. Ecol.102 (6), 1528–1536. doi: 10.1111/1365-2745.12300

  • CrossRef
  • Google Scholar

• 14

  GeraldiN. R.KleinS. G.AntonA.DuarteC. M. (2020). A framework for experimental scenarios of global change in marine systems using coral reefs as a case study. R. Soc. Open Sci.7 (1), 191118. doi: 10.1098/rsos.191118

  • CrossRef
  • Google Scholar

• 15

  GrennanA. K. (2007). The role of trehalose biosynthesis in plants. Plant Physiol.144 (1), 3–5. doi: 10.1104/pp.104.900223

  • CrossRef
  • Google Scholar

• 16

  HansonA. D.BeaudoinG. A.McCartyD. R.GregoryJ. F.III (2016). Does abiotic stress cause functional B vitamin deficiency in plants? Plant Physiol.172 (4), 2082–2097. doi: 10.1104/pp.16.01371

  • CrossRef
  • Google Scholar

• 17

  Heslop-HarrisonJ. S.SchwarzacherT.LiuQ. (2023). Polyploidy: its consequences and enabling role in plant diversification and evolution. Ann. Bot.131, 1–10. doi: 10.1093/aob/mcac132

  • CrossRef
  • Google Scholar

• 18

  HobdayA. J.OliverE. C. J.Sen GuptaA.BenthuysenJ. A.BurrowsM. T.DonatM. G.et al. (2018). Categorizing and naming marine heatwaves. Oceanography31, 162–173. doi: 10.5670/oceanog.2018.205

  • CrossRef
  • Google Scholar

• 19

  HoffmannA. A.SgròC. M. (2011). Climate change and evolutionary adaptation. Nature470 (7335), 479–485. doi: 10.1038/nature09670

  • CrossRef
  • Google Scholar

• 20

  HughesT. P.KerryJ. T.Álvarez-NoriegaM.Álvarez-RomeroJ. G.AndersonK. D.BairdA. H.et al. (2017). Global warming and recurrent mass bleaching of corals. Nature543 (7645), 373–377. doi: 10.1038/nature21707

  • CrossRef
  • Google Scholar

• 21

  IPCC. (2021). Summary For Policymakers. In: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Eds. Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S.L., Péan, C., Berger, S., Caud, N., Chen, Y., et al., pp. 3-32 Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press. 10.1017/9781009157896.001. Available at: https://www.ipcc.ch/report/ar6/wg1/.

  • Google Scholar

• 22

  JungE. M. U.Abdul MajeedN. A. B.BoothM. W.AustinR.SinclairE. A.FraserM. W.et al. (2023). Marine heatwave and reduced light scenarios cause species-specific metabolomic changes in seagrasses under ocean warming. New Phytol.239 (5), 1692–1706. doi: 10.1111/nph.19092

  • CrossRef
  • Google Scholar

• 23

  KendrickG. A.NowickiR. J.OlsenY. S.StrydomS.FraserM. W.SinclairE. A.et al. (2019). A systematic review of how multiple stressors from an extreme event drove ecosystem-wide loss of resilience in an iconic seagrass community. Front. Mar. Sci.6, 455. doi: 10.3389/fmars.2019.00455

  • CrossRef
  • Google Scholar

• 24

  KeuskampD. H.SasidharanR.VosI.PeetersA. J.VoesenekL. A.PierikR. (2011). Blue-light-mediated shade avoidance requires combined auxin and brassinosteroid action in Arabidopsis seedlings. Plant J.67 (2), 208–217. doi: 10.1111/j.1365-313X.2011.04597.x

  • CrossRef
  • Google Scholar

• 25

  KimD.PaggiJ. M.ParkC.BennettC.SalzbergS. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol.37 (8), 907–915. doi: 10.1038/s41587-019-0201-4

  • CrossRef
  • Google Scholar

• 26

  KingN. G.McKeownN. J.SmaleD. A.MooreP. J. (2018). The importance of phenotypic plasticity and local adaptation in driving intraspecific variability in thermal niches of marine macrophytes. Ecography41 (9), 1469–1484. doi: 10.1111/ecog.03186

  • CrossRef
  • Google Scholar

• 27

  KoskinenP.TörönenP.Nokso-KoivistoJ.HolmL. (2015). PANNZER: high-throughput functional annotation of uncharacterized proteins in an error-prone environment. Bioinformatics31 (10), 1544–1552. doi: 10.1093/bioinformatics/btu851

  • CrossRef
  • Google Scholar

• 28

  KrasenskyJ.JonakC. (2012). Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory networks. J. Exp. Bot.63 (4), 1593–1608. doi: 10.1093/jxb/err460

  • CrossRef
  • Google Scholar

• 29

  LaiX. J.YangR.LuoQ. J.ChenJ. J.ChenH. M.YanX. J. (2015). Glycerol-3-phosphate metabolism plays a role in stress response in the red alga Pyropia haitanensis. J. Phycol.51 (2), 321–331. doi: 10.1111/jpy.12276

  • CrossRef
  • Google Scholar

• 30

  LangfelderP.HorvathS. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinf.9 (1), 1–13. doi: 10.1186/1471-2105-9-559

  • CrossRef
  • Google Scholar

• 31

  LankauR.JørgensenP. S.HarrisD. J.SihA. (2011). Incorporating evolutionary principles into environmental management and policy. Evol. Appl.4 (2), 315–325. doi: 10.1111/j.1752-4571.2010.00171.x

  • CrossRef
  • Google Scholar

• 32

  LefcheckJ. S.WilcoxD. J.MurphyR. R.MarionS. R.OrthR. J. (2017). Multiple stressors threaten the imperiled coastal foundation species eelgrass (Zostera marina) in Chesapeake Bay, USA. Global Change Biol.23 (9), 3474–3483. doi: 10.1111/gcb.13623

  • CrossRef
  • Google Scholar

• 33

  LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics25 (16), 2078–2079. doi: 10.1093/bioinformatics/btp352

  • CrossRef
  • Google Scholar

• 34

  LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15 (12), 1–21. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar

• 35

  MalandrakisE.DadaliO.KavourasM.DanisT.PanagiotakiP.MiliouH.et al. (2017). Identification of the abiotic stress-related transcription in little Neptune grass Cymodocea nodosa with RNA-seq. Mar. Genomics34, 47–56. doi: 10.1016/j.margen.2017.03.005

  • CrossRef
  • Google Scholar

• 36

  Marín-GuiraoL.RuizJ. M.DattoloE.Garcia-MunozR.ProcacciniG. (2016). Physiological and molecular evidence of differential short-term heat tolerance in Mediterranean seagrasses. Sci. Rep.6 (1), 28615. doi: 10.1038/srep28615

  • CrossRef
  • Google Scholar

• 37

  MoriyaY.ItohM.OkudaS.YoshizawaA. C.KanehisaM. (2007). KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res.35, W182–W185. doi: 10.1093/nar/gkm321

  • CrossRef
  • Google Scholar

• 38

  NitschkeS.CortlevenA.IvenT.FeussnerI.HavauxM.RieflerM.et al. (2016). Circadian stress regimes affect the circadian clock and cause jasmonic acid-dependent cell death in cytokinin-deficient Arabidopsis plants. Plant Cell28 (7), 1616–1639. doi: 10.1105/tpc.16.00016

  • CrossRef
  • Google Scholar

• 39

  OliverE. C.DonatM. G.BurrowsM. T.MooreP. J.SmaleD. A.AlexanderL. V.et al. (2018). Longer and more frequent marine heatwaves over the past century. Nat. Commun.9 (1), 1–12. doi: 10.1038/s41467-018-03732-9

  • CrossRef
  • Google Scholar

• 40

  OomenR. A.HutchingsJ. A. (2017). Transcriptomic responses to environmental change in fishes: insights from RNA sequencing. Facets2 (2), 610–641. doi: 10.1139/facets-2017-0015

  • CrossRef
  • Google Scholar

• 41

  Pantoja-BenavidesA. D.Garces-VaronG.Restrepo-DíazH. (2022). Foliar cytokinins or brassinosteroids applications influence the rice plant acclimatization to combined heat stress. Front. Plant Sci.13, 983276. doi: 10.3389/fpls.2022.983276

  • CrossRef
  • Google Scholar

• 42

  PeclG. T.AraújoM. B.BellJ. D.BlanchardJ.BonebrakeT. C.ChenI. C.et al. (2017). Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. Science355 (6332), eaai9214. doi: 10.1126/science.aai9214

  • CrossRef
  • Google Scholar

• 43

  PerteaG.PerteaM. (2020). GFF utilities: GffRead and GffCompare. F1000Research9, 304. doi: 10.12688/f1000research.23297.1

  • CrossRef
  • Google Scholar

• 44

  PerteaM.PerteaG. M.AntonescuC. M.ChangT. C.MendellJ. T.SalzbergS. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol.33 (3), 290–295. doi: 10.1038/nbt.3122

  • CrossRef
  • Google Scholar

• 45

  RendeU.NiittyläT.MoritzT. (2019). Two-step derivatization for determination of sugar phosphates in plants by combined reversed phase chromatography/tandem mass spectrometry. Plant Methods15 (1), 1–10. doi: 10.1186/s13007-019-0514-9

  • CrossRef
  • Google Scholar

• 46

  RiveraH. E.AichelmanH. E.FiferJ. E.KriefallN. G.WuitchikD. M.SmithS. J.et al. (2021). A framework for understanding gene expression plasticity and its influence on stress tolerance. Mol. Ecol.30 (6), 1381–1397. doi: 10.1111/mec.15820

  • CrossRef
  • Google Scholar

• 47

  RoeberV. M.BajajI.RohdeM.SchmüllingT.CortlevenA. (2021). Light acts as a stressor and influences abiotic and biotic stress responses in plants. Plant Cell Environ.44 (3), 645–664. doi: 10.1111/pce.13948

  • CrossRef
  • Google Scholar

• 48

  SadakM. S.El-BassiounyH. M. S.DawoodM. G. (2019). Role of trehalose on antioxidant defense system and some osmolytes of quinoa plants under water deficit. Bull. Natl. Res. Centre43 (1), 1–11. doi: 10.1186/s42269-018-0039-9

  • CrossRef
  • Google Scholar

• 49

  SainiS.SharmaI.PatiP. K. (2015). Versatile roles of brassinosteroid in plants in the context of its homoeostasis, signaling and crosstalks. Front. Plant Sci.6, 950. doi: 10.3389/fpls.2015.00950

  • CrossRef
  • Google Scholar

• 50

  Santillán-SarmientoA.PazzagliaJ.RuoccoM.DattoloE.AmbrosinoL.WintersG.et al. (2023). Gene co-expression network analysis for the selection of candidate early warning indicators of heat and nutrient stress in Posidonia oceanica. Sci. Total Environ.877, 162517. doi: 10.1016/j.scitotenv.2023.162517

  • CrossRef
  • Google Scholar

• 51

  ScheffersB. R.De MeesterL.BridgeT. C.HoffmannA. A.PandolfiJ. M.CorlettR. T.et al. (2016). The broad footprint of climate change from genes to biomes to people. Science354 (6313), aaf7671. doi: 10.1126/science.aaf7671

  • CrossRef
  • Google Scholar

• 52

  SewelamN.KazanK.SchenkP. M. (2016). Global plant stress signaling: reactive oxygen species at the cross-road. Front. Plant Sci.7, 187. doi: 10.3389/fpls.2016.00187

  • CrossRef
  • Google Scholar

• 53

  SinclairE. A.HoveyR. K.KraussS. L.AnthonyJ. M.WaycottM.KendrickG. A. (2023). Historic and contemporary biogeographic perspectives on range-wide spatial genetic structure in a widespread seagrass. Ecol. Evol.13, e9900. doi: 10.1002/ece3.9900

  • CrossRef
  • Google Scholar

• 54

  SinghA. P.Savaldi-GoldsteinS. (2015). Growth control: brassinosteroid activity gets context. J. Exp. Bot.66 (4), 1123–1132. doi: 10.1093/jxb/erv026

  • CrossRef
  • Google Scholar

• 55

  Sobieszczuk-NowickaE. (2017). Polyamine catabolism adds fuel to leaf senescence. Amino Acids49 (1), 49–56. doi: 10.1007/s00726-016-2377-y

  • CrossRef
  • Google Scholar

• 56

  StillmanJ. H.ArmstrongE. (2015). Genomics are transforming our understanding of responses to climate change. BioScience65 (3), 237–246. doi: 10.1093/biosci/biu219

  • CrossRef
  • Google Scholar

• 57

  StrydomS.MurrayK.WilsonS.HuntleyB.RuleM.HeithausM.et al. (2020). Too hot to handle: Unprecedented seagrass death driven by marine heatwave in a World Heritage Area. Global Change Biol.26 (6), 3525–3538. doi: 10.1111/gcb.15065

  • CrossRef
  • Google Scholar

• 58

  SupekF.BošnjakM.ŠkuncaN.ŠmucT. (2011). REVIGO summarizes and visualizes long lists of gene ontology terms. PloS One6 (7), e21800. doi: 10.1371/journal.pone.0021800

  • CrossRef
  • Google Scholar

• 59

  TohgeT.WatanabeM.HoefgenR.FernieA. R. (2013). Shikimate and phenylalanine biosynthesis in the green lineage. Front. Plant Sci.4, 62. doi: 10.3389/fpls.2013.00062

  • CrossRef
  • Google Scholar

• 60

  TsioliS.KoutalianouM.GkafasG. A.ExadactylosA.PapathanasiouV.KatsarosC. I.et al. (2022). Responses of the Mediterranean seagrass Cymodocea nodosa to combined temperature and salinity stress at the ionomic, transcriptomic, ultrastructural and photosynthetic levels. Mar. Environ. Res.175, 105512. doi: 10.1016/j.marenvres.2021.105512

  • CrossRef
  • Google Scholar

• 61

  Van de PeerY.AshmanT.-L.SoltisP. S.SoltisD. E. (2021). Polyploidy: an evolutionary and ecological force in stressful times. Plant Cell33 (1), 11–26. doi: 10.1093/plcell/koaa015

  • CrossRef
  • Google Scholar

• 62

  WernbergT.BennettS.BabcockR. C.De BettigniesT.CureK.DepczynskiM.et al. (2016). Climate-driven regime shift of a temperate marine ecosystem. Science353 (6295), 169–172. doi: 10.1126/science.aad8745

  • CrossRef
  • Google Scholar

• 63

  WickhamH.WickhamH. (2016). Data Analysis (Cham, Switzerland: Springer International Publishing), 189–201.

  • Google Scholar

• 64

  WijffelsS. E.BeggsH.GriffinC.MiddletonJ. F.CahillM.KingE.et al. (2018). A fine spatial-scale sea surface temperature atlas of the Australian regional seas (SSTAARS): Seasonal variability and trends around Australasia and New Zealand revisited. J. Mar. Syst.187, 156–196. doi: 10.1016/j.jmarsys.2018.07.005

  • CrossRef
  • Google Scholar

• 65

  WuT.HuE.XuS.ChenM.GuoP.DaiZ.et al. (2021). clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation2 (3), 100141. doi: 10.1016/j.xinn.2021.100141

  • CrossRef
  • Google Scholar

• 66

  YooH.WidhalmJ. R.QianY.MaedaH.CooperB. R.JannaschA. S.et al. (2013). An alternative pathway contributes to phenylalanine biosynthesis in plants via a cytosolic tyrosine: phenylpyruvate aminotransferase. Nat. Commun.4 (1), 2833. doi: 10.1038/ncomms3833

  • CrossRef
  • Google Scholar

• 67

  YorkP. H.GruberR. K.HillR.RalphP. J.BoothD. J.MacreadieP. I. (2013). Physiological and morphological responses of the temperate seagrass Zostera muelleri to multiple stressors: investigating the interactive effects of light and temperature. PloS One8 (10), e76377. doi: 10.1371/journal.pone.0076377

  • CrossRef
  • Google Scholar

• 68

  ZandalinasS. I.BalfagónD.Gómez-CadenasA.MittlerR. (2022). Plant responses to climate change: metabolic changes under combined abiotic stresses. J. Exp. Bot.73 (11), 3339–3354. doi: 10.1093/jxb/erac073

  • CrossRef
  • Google Scholar

• 69

  ZandalinasS. I.MittlerR. (2022). Plant responses to multifactorial stress combination. New Phytol.234 (4), 1161–1167. doi: 10.1111/nph.18087

  • CrossRef
  • Google Scholar

• 70

  ZhangH.SonnewaldU. (2017). Differences and commonalities of plant responses to single and combined stresses. Plant J.90 (5), 839–855. doi: 10.1111/tpj.13557

  • CrossRef
  • Google Scholar

• 71

  ZhangJ.YangQ.YueW.YangB.ZhouW.ChenL.et al. (2023). Seagrass Thalassia hemprichii and associated bacteria co-response to the synergistic stress of ocean warming and ocean acidification. Environ. Res.236, 116658. doi: 10.1016/j.envres.2023.116658

  • CrossRef
  • Google Scholar

• 72

  ZhangY.ZhaoP.YueS.LiuM.QiaoY.XuS.et al. (2021). New insights into physiological effects of anoxia under darkness on the iconic seagrass Zostera marina based on a combined analysis of transcriptomics and metabolomics. Sci. Total Environ.768, 144717. doi: 10.1016/j.scitotenv.2020.144717

  • CrossRef
  • Google Scholar