Wild G. bailinae specimens were collected from Wushi (N 20°33′, E 109°51′), Guangdong Province, China, in July 2021.The samples were brought back to the lab in a portable fridge at 4°C. Fresh algae were randomly selected and cleaned thoroughly using sterilized seawater. They were then temporarily cultivated in a tank filled with sterile seawater that had been enhanced with f/2 media at 25°C under a 12 L:12 D photoperiod and light intensity of 100 μmol·m⁻²·s⁻¹. Three beakers containing 1 L of fresh sterile seawater enriched with f/2 medium were prepared and 3 g of healthy algae was added to each of them. The beakers were then transferred to incubators set at 15°C (LT: low temperature), 25°C (MT: middle temperature), and 35°C (HT: high temperature) for 7 days. Three biological replicates were prepared per treatment and the medium was changed every 2 days. After 7 days of culture, the algae were used for growth and F _v/F _m measurements, as well as transcriptomic analysis.

After the treatment, the fresh weight of G. bailinae was measured and the relative growth rate (RGR) was determined using the formula supplied by Yong (Yong et al., 2013):

where M ₀ represents the starting weight of the algae and M ₇ represents their weight after the 7-day treatment.

The maximum photochemical quantum yield of PSII (F _v/F _m) was measured based on Li et al. (Li et al., 2016) after 7 days of treatment at different temperatures using a Hansatech FMS-2 modulated fluorometer (Hansatech Instruments, Norfolk, UK). The formula used was F _v = F _m – F ₀; the minimum (F ₀) and maximum (F _m) fluorescence were obtained by exposing algae that had adapted to darkness for 30 minutes to actinic light and saturated light provided by the fluorometer (Li et al., 2016).

Total RNA was isolated from G. bailinae using the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) in accordance with the manufacturer’s instructions after 7 days of LT, MT, and HT treatments. An Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) was used to analyze the quality of total RNA and the results were validated via RNase-free agarose gel electrophoresis. Following total RNA extraction, Oligo (dT) beads were used to enrich G. bailinae’s mRNA. Using random primers, the enriched mRNA was reverse transcribed into cDNA after being broken into small fragments by applying fragmentation buffer. Subsequently, DNA polymerase I, RNase H, dNTP, and buffer were used to create second-strand cDNA. Following purification with Qia Quick PCR extraction kit (Qiagen, Venlo, Netherlands), the cDNA fragments were ligated to the Illumina sequencing adapter after end repair and poly(A) addition. The ligation products were size selected via agarose gel electrophoresis, amplified via PCR, and then sequenced using Illumina NovaSeq 6000 (Gene Denovo Biotechnology Co., Guangzhou, China).

Principal component analysis (PCA) was carried out using R (http://www.r-project.org/) on all nine transcriptome datasets to investigate variations in expression patterns between samples. The three duplicates of each temperature treatment were subjected to Pearson correlation analysis in order to evaluate the reproducibility of the transcriptome data.

The raw reads were trimmed by removing those containing adapters, more than 10% of unknown nucleotides, and more than 50% of low quality (Q-value ≤ 20) bases using fastp (version 0.18.0) (Chen et al., 2018). The transcriptome’s de novo assembly was carried out using the short-read assembling program Trinity (Grabherr et al., 2011). BLASTx (http://www.ncbi.nlm.nih.gov/BLAST/) with an E-value threshold of 10⁵ was used to compare transcripts to the NCBI non-redundant protein (Nr) database (http://www.ncbi.nlm.nih.gov), the Swiss-Prot protein database (http://www.expasy.ch/sprot), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg), and the COG/KOG database (http://www.ncbi.nlm.nih.gov/COG), and annotate the unigenes. Protein function annotations were then obtained based on the best alignment results.

The differential expression analysis of RNAs between different groups and between samples was performed using DESeq2 (Love et al., 2014) and edgeR (Robinson et al., 2010) software, respectively. The genes with a parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≥ 3 were considered DEGs.

Nine genes associated with temperature changes were randomly selected from the transcriptome results to validate the transcriptome using the 18sRNA (Liu et al., 2021a) gene as a reference gene. The following reaction steps were adopted: 95°C for 300 s, followed by 45 cycles at 95°C for 10 s, 60°C for 15 s, 72°C for 15 s, and reading the fluorescence signal, followed by 1 cycle at 95°C for 10 s, 65°C for 60 s, and 95°C for 1 s. Each reaction was conducted in triplicate. Primer Premier 5.0 (Premier, Toronto, ON, Canada) was used to design each primer. The 2^−ΔΔCt (RQ) method was used to determine the relative transcript level (Qin et al., 2021). Primer sequences and qRT-PCR validation results are included in Supplementary Table S1 .

The RGR and F _v/F _m data were expressed as means ± SD (n ≥ 3). The statistical significance of different treatments was assessed using one-way ANOVA and Duncan’s test in SPSS 25 (IBM Corp., Armonk, NY, USA) with a significance level of P< 0.05.

G. bailinae maintained a positive growth during the period of culture under the three different temperatures. The HT treatment (35°C) resulted in a maximum growth rate of 2.1%, which was significantly higher than the values obtained at 15°C and 25°C (P< 0.01); in particular, it was 15 times higher than that obtained at 15°C (0.14%) ( Figure 1A ). The F _v/F _m values at 25°C and 35°C were 0.65 and 0.62, respectively, which were significantly higher than the value at 15°C (0.36) (P< 0.01) ( Figure 1B ).

The raw data obtained from the nine samples sequenced were subjected to quality control to obtain clean data ( Supplementary Table S2 ). Out of a total of 494,707,094 reads 492,443,386 were retained after removing adapters and low-quality sequences. In each sample, the yield of clean reads was above 99%, the GC content was between 50.94% and 51.56%, and the Q20(%) value was above 98%. The clustering of the assembled transcripts of G. bailinae yielded a total of 20,049 single genes (> 200 base pairs) with an N50 count of 4,140 and a length of 1,072 base pairs. The size of single genes ranged from 201 to 8,602 base pairs, with an average length of 769 base pairs and a total length of 15,419,753 base pairs ( Supplementary Table S3 ).

The three duplicate samples in each treatment had Pearson correlation coefficients (r) of less than 0.96 (r > 0.80 indicates a strong correlation) ( Supplementary Figure S1A ). PCA showed that the expression patterns of G. bailinae differed under the three temperatures ( Supplementary Figure S1B ). Furthermore, qRT-PCR and mRNA-Seq analyses showed that nine randomly selected genes had a similar expression tendency ( Supplementary Table S1 ).

At present, the whole-genome sequencing of G. bailinae has not been completed, therefore, only 20,049 gene transcript sequences were compared to the NCBI non-redundant protein databases Nr, Swiss-Prot, KEGG, and COG/KOG (evalue< 0.00001) using BLASTx. The proteins with the highest sequence similarity to the given unigene and its functional annotation information were obtained. In total, 20,049 genes were assembled from the transcriptome data in this study and 13,095 (65.31%) of them were compared to the four major databases. Of these 20,049 genes, 12,828 (63.98%), 5,246 (26.17%), 5,027 (25.07%), and 9,714 (48.45%) genes matched the Nr, KEGG, KOG, and GO databases, respectively.

In the Nr database, the main match obtained was Gracilariopsis chorda (84.82%), which belongs to the same genus as G. bailinae, in line with the taxonomic status distinction ( Supplementary Figure S2 ). Cellular processes, environmental information processing, genetic information processing, metabolism, and organic systems were the five primary functional categories of the genes that matched the KEGG database. Each category was subdivided into secondary groups of functions. In particular, the largest number of unigenes were involved in “global and overview maps” (1,445), followed by “translation” (750) and “folding, sorting, and degradation” (409) ( Supplementary Figure S3 ). The genes that matched the KOG database were divided into 25 groups: the largest was “translation, ribosomal structure and biogenesis” (781), followed by “general function prediction only” (701), and the smallest was “cell motility” (5) ( Supplementary Figure S4 ). A total of 5,824 (29.05%) genes were matched to Swiss-Prot database. The number of genes matched in all databases was 4611 (23.00%) ( Supplementary Figure S5 ).

Based on the Nr annotation information, Blast2GO was used to obtain the GO function annotations (Ashburner et al., 2000; Conesa et al., 2005). Three ontologies in GO were used to categorize the genes’ molecular functions (MF), cellular components (CC), and biological processes (BP). The 9,714 unigenes annotated into the GO database by GO annotation were classified into 50 subclasses, accounting for 48.45% of all unigenes, and the same unigenes could be annotated to different subclasses. The top three subclasses of the GO ontologies were: cellular process, metabolic process, and localization with 6,637, 5,845, and 1,497, genes, respectively, for BP; cellular anatomical entity, protein-containing complex, and virion component, with 4,253, 1,958 and 110 genes, respectively, for CC; and catalytic activity, binding, and transporter activity, with 5,066, 4,771, and 802 genes, respectively, for MF ( Supplementary Figure S6 ).

Macroalgae are expected to be especially sensitive to ocean warming since the rates of their biochemical and physiological processes are substantially regulated by ambient temperature (as observed in ectotherms) (Davison, 1991; Brown et al., 2004; Graba-Landry et al., 2018). Typically, within a species’ thermal tolerance range, the rates of physiological processes and other parameters (e.g., growth) increase with temperature until they reach a thermal optimum and rapidly decline thereafter (Brown et al., 2004; Graba-Landry et al., 2018). Different algae show different adaptations to temperature variation, and growth rate is an important indicator that visually reflects their growth status (Yang et al., 2021; Gao et al., 2022): the highest growth rate is recorded at the optimum growth temperature. For example, the highest growth rate of G. lemaneiformis (Liu et al., 2017), Ulva prolifera (Cui et al., 2015), and Caulerpa sertularioides (Zhong et al., 2021) was observed at 23°C, 20°C, and 26°C, respectively. In the present study, the RGR of G. bailinae reached its maximum at 35°C and growth was significantly inhibited under the lower temperatures tested. F _v/F _m is the maximum photochemical quantum yield of PSII, which reflects the potential maximum photosynthetic capacity of macroalgae and can also indirectly reflect their growth state and living environment (Goltsev et al., 2016; Cui et al., 2019; Zhong et al., 2021). Studies have indicated that F _v/F _m is a heat-sensitive indicator, which makes it an important tool for the evaluation of plants’ heat tolerance (Yamada et al., 1996; Wahid et al., 2007; Wu et al., 2019). For example, in G. lemaneiformis, this parameter was shown to decrease under unsuitable temperatures (Liu et al., 2017). A notable drop in F _v/F _m was also observed in Saccharina japonica under high temperature stress, which was caused by the deactivation of PSII (Shindo et al., 2022). In the present study, the analysis of F _v/F _m in G. bailinae indicated that 25°C and 35°C were suitable temperatures, but the parameter decreased significantly at 15°C. Obviously, the changes in RGR and F _v/F _m under different temperatures indicate that G. bailinae has a remarkable capacity to adapt to high temperatures.

It has long been acknowledged that plant photosynthesis is one of the most temperature-sensitive processes (Li et al., 2019). This process depends on electron transport into the photosynthetic system and it gives algae the energy they need to grow and develop (Qin et al., 2021). In the cells of red algae, light energy absorbed by the phycobillisome and the light-harvesting chlorophyll protein complex (LHC) drives the flow of electrons from H₂O to NADP⁺, which promotes the reduction of NADP⁺ to NADPH (Croce and Van Amerongen, 2014). At the same time, ATP synthase uses the electrochemical gradient generated simultaneously across the thylakoid membrane for ATP synthesis (Croce and Van Amerongen, 2014). NADPH and ATP represent the basis for all subsequent chemical reactions (Croce and Van Amerongen, 2014). Algae can regulate the above processes by controlling gene expression to resist or adapt to environmental changes, that is, regulating the synthesis of NADPH and ATP. Wang et al. (2018) reported that the protein-related genes of the photosynthetic antenna in wild-type P. haitanensis were significantly downregulated under high temperature stress, impeding the absorption of sufficient light to activate photosynthesis, and therefore resulting in insufficient energy to resist heat stress (Wang et al., 2018). Another study on Amphiroa fragilissima revealed that the expression levels of several DEGs involved in photosynthesis were suppressed under high temperature stress compared to those in the control group (Yang et al., 2021). However, in the present study, the photosynthetic activity of G. bailinae remained normal despite the downregulation of photosynthesis-related genes, which was not in line with the results of the above-mentioned studies. The downregulation of genes in most heat-sensitive algae under temperature stress may lead to fatal results, but G. bailinae has a high temperature adaptability, and 35°C is its suitable growth temperature. A study on the high-temperature tolerant P. haitanensis showed that this species can reduce unnecessary energy consumption by decreasing the rate of photosynthesis and energy metabolism (Wenlei et al., 2018). A proteomic analysis similarly indicated that the P. haitanensis strain Z-61 resists high temperature stress by downregulating photosynthesis (Xu et al., 2014). Similarly, G. bailinae regulates NADPH and ATP synthesis by downregulating photosynthesis-related genes to decrease the electron transfer rate and H⁺ production in the photosynthetic system, thereby potentially reducing unnecessary energy consumption and ultimately adapting to high temperature environments. In addition, as observed for the compensatory reaction of gene expression (Dal Bosco et al., 2004), the high expression of LHCA1-related genes in photosynthetic antenna proteins may compensate for the decreased light capture efficiency caused by the downregulation of genes related to LHCB2, apcB, and cpcE, so as to meet the light energy demand of G. bailinae and achieve a higher growth rate.

In plants, a metabolic regulatory network is frequently engaged in the response to abiotic stimuli (Wenlei et al., 2018). Previous studies have shown that P. haitanensis under stress can reduce unnecessary energy consumption to prevent cell damage (Xu et al., 2014; Xu et al., 2016). However, when stress levels become excessive, the corresponding energy shortage in algae would result in an enhancement of the inherent pathways of carbohydrate metabolism and the induction of alternative pathways, such as glycolysis, to provide energy and the related carbon skeletons for key metabolic processes (Obata and Fernie, 2012; Wenlei et al., 2018). Wang et al. (Wang et al., 2018) corroborated this claim by studying the heat-tolerant P. haitanensis, whose genes showed a tendency to be upregulated under high temperature stress to maintain the energy and associated carbon skeleton for key metabolic processes. In the present study, a joint analysis was performed for the following five pathways related to the carbohydrate and energy metabolisms: Calvin cycle, glycolysis/gluconeogenesis, PPP, pyruvate metabolism, and TCA cycle. In our study, we found that only the DEG encoding GPI was significantly upregulated under HT treatment, while other DEGs enriched in these pathways were downregulated ( Figure 4 ). The accumulation of sugar phosphates associated with glycolysis and the oxidative pentose phosphate pathway (OPPP) indicated that glycolytic carbon flow may be rerouted into the OPPP to produce NADPH for antioxidative action (Obata and Fernie, 2012; Wenlei et al., 2018). GPI catalyzes the interconversion between F6P and G6P. Further oxidation of G6P in the PPP can produce a large amount of NADPH, which may be one of the remedies for the decrease in NADPH synthesis during photosynthesis. The Calvin cycle, which comprises 13 chemical steps and 11 enzymes that catalyze them, is the main mechanism for carbon fixation in C3 plants (Tamoi et al., 2005). The cycle converts atmospheric CO₂ into carbon skeletons that are utilized to produce starch and sucrose using the byproducts of the light reactions of photosynthesis (i.e., ATP and NADPH) (Raines, 2003). Gene downregulation during the Calvin cycle is generally thought to lead to a reduced synthesis of related enzyme proteins, especially the key protein Rubisco. However, in this study, the expression of genes related to the Rubisco small subunit rbcS was downregulated under HT treatment, but the biomass of G. bailinae was relatively high. Related research shows that Rubisco activase induced the activation of Rubisco, but the synthesis of this protein (transcription and translation) and enzyme had different sensitivities to long-term heat stress (Almeselmani et al., 2012). Therefore, due to the adaptability of G. bailinae to high temperatures, Rubisco activase may activate Rubisco more obviously under high temperature conditions, thereby improving the carbon sequestration efficiency of the protein and reducing its absolute content during the Calvin cycle. However, further studies are needed to confirm this. The TCA cycle promotes the oxidation of respiratory substrates for ATP generation, and therefore plays a crucial role in plant resistance to abiotic challenges (Sweetlove et al., 2010; Obata and Fernie, 2012; Wenlei et al., 2018). It is worth noting that the genes encoding ACLY and IDH3 were significantly downregulated in the TCA cycle. IDH3 catalyzes the conversion of isocitrate to alpha-ketoglutarate to produce NADH. The downregulation of IDH3 indicated that G. bailinae adapts to high temperatures by reducing unnecessary energy consumption. In addition, the downregulation of other genes in the Calvin cycle, glycolysis/gluconeogenesis, and pyruvate metabolism pathways also indicated that G. bailinae slows the rate of energy or carbohydrate metabolisms to adapt to high temperatures. On the contrary, under low temperature stress, due to energy shortage, this macroalgae needs to produce more energy by overexpressing genes associated with the synthesis of energy molecules, such as MDH1 and SDHB, to maintain the key life activities of cells.

In response to temperature stress, plants typically control the percentage of polyunsaturated fatty acids to alter membrane fluidity (Sun et al., 2015). Two aspects are involved in the production of unsaturated fatty acids: fatty acid peptide chain extension and desaturation (Qin et al., 2021). Fatty acids, which are the major components of membrane phospholipids, induce desaturation to promote cell membrane fluidity and stability, hence minimizing cellular damage (Qin et al., 2021). Transcriptome analysis of G. lemaneiformis revealed that the expression levels of genes associated with fatty acid extension and desaturation under normal temperatures were significantly lower than those under low temperature stress (Qin et al., 2021). Wang et al. (2018) also pointed out that genes related to unsaturated fatty acid synthesis in P. haitanensis were more significantly downregulated under normal temperatures than under high temperature stress. In this study, the genes associated with desaturation (SCD) were significantly downregulated under the HT treatment, while those associated with both fatty acid chain extension (KAR) and desaturation (SCD) were significantly upregulated under the LT treatment. This indicates that G. bailinae may regulate the membrane’s lipid fluidity to adapt to high temperatures by gradually downregulating the genes related to the biosynthesis of unsaturated fatty acids.

Alpha-linolenic acid (a C18 polyunsaturated acid) is a biosynthetic substrate for jasmonic acid (JA) and its diverse derivatives, which are lipid-derived signaling molecules involved in the regulation of numerous plant processes, including growth, reproductive development, photosynthesis, and responses to abiotic and biotic stresses (Bhatla, 2018). In plants, alpha-linolenic acid is perhaps the most prevalent and chemically diverse source of polyunsaturated fatty acids (Gerwick et al., 1991; Leung et al., 2022). JA can further form methyl jasmonate, which is considered to be crucial for algae’s ability to withstand heat (Qin et al., 2021). Exogenous methyl jasmonate has been found to reduce the negative effects of high temperature stress on a number of pathways in G. lemaneiformis (Wang et al., 2017). A transcriptome analysis of this species under low temperature stress also revealed that the activation of crucial genes associated with methyl jasmonate production was dramatically elevated (Qin et al., 2021). Su et al. (Su et al., 2019) reported that heat stress induced the upregulation of genes involved in the metabolism of alpha-linolenic acid in wheat. In the present study, it was revealed that the genes encoding OPR (the enzyme 12-oxophytodienoic acid reductase that catalyzes the conversion of 12-OPDA to OPC8) were significantly downregulated under the HT treatment. In contrast, those encoding OPR and ACX (acyl-CoA oxidase, which plays a crucial role in the three rounds of β-oxidation for the synthesis of JA acid from alpha-linolenic acid (Liu et al., 2021b)) were significantly upregulated in the LT treatment ( Figure 6 ). These results suggest that G. bailinae can adapt to high temperatures by gradually limiting JA synthesis.