Transcriptome analysis of Sonneratia caseolaris seedlings under chilling stress

Sonneratia caseolaris is a native mangrove species found in China. It is fast growing and highly adaptable for mangrove afforestation, but suffered great damage by chilling event once introduced to high latitude area. To understand the response mechanisms under chilling stress, physiological and transcriptomic analyses were conducted. The relative electrolyte conductivity, malondialdehyde (MDA) content, soluble sugar content and soluble protein content increased significantly under chilling stress. This indicated that S. caseolaris suffered great damage and increased the levels of osmoprotectants in response to the chilling stress. Gene expression comparison analysis of S. caseolaris leaves after 6 h of chilling stress was performed at the transcriptional scale using RNA-Seq. A total of 168,473 unigenes and 3,706 differentially expressed genes (DEGs) were identified. GO and KEGG enrichment analyses showed that the DEGs were mainly involved in carbohydrate metabolism, antioxidant enzyme, plant hormone signal transduction, and transcription factors (TFs). Sixteen genes associated with carbohydrate metabolism, antioxidant enzyme, phytohormones and TFs were selected for qRT-PCR verification, and they indicated that the transcriptome data were reliable. Our work provided a comprehensive review of the chilling response of S. caseolaris at both physiological and transcriptomic levels, which will prove useful for further studies on stress-responses in mangrove plants.


INTRODUCTION
The growth, development, and geographical distribution of plants are limited by temperatures in their habitat (Jin & Kim, 2013).Low temperature stress can be divided into chilling stress (<10 • C) and freezing stress (<0 • C) (Zhou et al., 2011).Plants have evolved a series of response mechanisms to adapt to low-temperature environments.These include regulation of gene expression and physiological and biochemical process changes that enhance cold resistance (Thomashow, 1999;Sung et al., 2003), These changes have been studied in model plants like Arabidopsis thaliana and some cereal crops (Buti et al., 2018;Zhang et al., 2004;Román et al., 2012).Response mechanisms enhancing cold resistance To understand the mechanisms of S. caseolaris under chilling stress, we first measured the physiological changes in S. caseolaris under chilling stress, Then, we used RNA-seq for comparative analyses of S. caseolaris leaves between a control group (0 h, CK) and a chilling treated group (4 • C 6 h, CT).The reliability of the transcriptome results was verified by real-time quantitative polymerase chain reaction (qRT-PCR).The data improved our understanding of the chilling response mechanism in S. caseolaris and identified potential cold tolerance gene targets in this, and other mangrove species.

Plant materials and treatments
We collected ripe fruits of S. caseolaris in Hainan Dongzhaigang National Nature Reserve (19 • 57 23 N, 110 • 30 10 E) in China.The fruit was soaked in distilled water for 7 d and then individually planted in pots with a mud and sand (1:1) mixture as support.Growth chamber conditions were 30/28 • C (day/night) and a 14:10 h (L: D) photoperiod, 400 µ mol m −2 s −1 photon flux density, and 90% constant humidity.The germinated seeds were cultured in a growth chamber for 5 weeks, at which time the seedlings were about 12 cm high with 12 leaves.We performed the chilling stress treatment on these seedlings in a growth chamber, by setting the temperature to 4 • C and leaving the other conditions unchanged.The second leaf from the top of the seedlings were collected as samples from 0 h (CK) 3 h, 6 h, 12 h to 24 h separately, and then frozen in liquid nitrogen and stored at −80 • C. Each treatment had three biological replicates.

Physiological response assays of chilling-treated seedlings
We measured relative electrolyte conductivity following Wang et al. (2017).The malondialdehyde (MDA) content was measured by the method of Yao et al. (2011).The soluble sugar content was measured by the method of Wang et al. (2013).The soluble protein content was measured by using the Bradford method (Bradford, 1976).Three biological replicates were measured at each time point.

RNA extraction
Based on the physiological results, we used the seedlings without chilling as the control group (CK_1, CK_2, CK_3) and seedlings under chilling stress for 6 h were set as the treatment group (CT_1, CT_2, CT_3) for RNA-seq.Total RNA of leaf samples was extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to manufacturer instructions.RNA purity was checked by T6 Spectrophotometer (Puxi, Beijing, China).RNA integrity and concentration were measured by the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

cDNA library construction and transcriptome sequencing
A 3 µg amount of total RNA of each sample was enriched using magnetic beads with Oligo (dT) and then fragmented into short fragments ∼200 bp by an interrupting reagent.The first strand cDNA was synthesized by a random hexamer primer, and the second strand of cDNAs was synthesized using the buffer, dNTPs, DNA polymerase I, and RNase H.After the double-stranded cDNA was generated and processed with end-repair, A-tailed, and adapters ligation, we constructed a cDNA library by PCR amplification.The cDNA library was then sequenced using the BGISEQ-500 platform with paired-end reads of 100 bp, by the Beijing Genomic Institution (www.genomics.org.cn, BGI, Shenzhen, China).All raw reads data from this study were submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive database under the accession number PRJNA668053.

Transcripts assembly and functional gene annotation
The raw data were filtered with SOAPnuke (v1.5.2) (Li et al., 2008).We discarded the low quality reads such as those with adapter sequences, ambiguous bases 'N' is more than 5%, with more than 20% Q <10 bases.The remaining clean data were used for de novo assembly with the Trinity program (v2.0.6) (Grabherr et al., 2011) and then the unigene sequences were aligned with the following databases: NR, NT, SwissProt, KEGG, KOG, Pfam and GO by BLAST software (version 2.2.26).

Differential expression analyses
The expression level of each transcript was measured by RSEM (v1.2.12), and the expression level of a single gene was expressed as the value of FPKM (fragments per kilobase of transcript per million mapped reads) method (Li & Dewey, 2011).For the sake of determining the differentially expressed genes (DEGs), we performed differential expression analysis by the DESeq2 (Anders & Huber, 2010).A P-value ≤ 0.05 and an expression level |log2(ratio)|>2 were regarded as significant differential expression.Gene ontology (GO) enrichment analysis of DEGs was conducted using the GOseq R package, GO terms with p < 0.05 were regarded as significantly enriched in DEGs (Young et al., 2010).We analyzed the KEGG pathway functional enrichment of DEGs using the R package.Significant enrichment was defined as when the FDR terms were < 0.001.

Quantitative real-time PCR analysis
We selected 16 DEGs to verify the reliability of transcriptome sequencing through qRT-PCR.The specific primers of qRT-PCR were designed by Primer3 software (http://primer3.ut.ee/) and are listed in Table S1.The qRT-PCR analyses were conducted using a LightCycler 480 RT-PCR system (Roche Applied Science, Mannheim, Germany) under the following conditions: 95 • C for 3 min, 40 cycles of 95 • C for 15 s, 56 • C for 30 s, and 72 • C for 20 s.Three replicates of each cDNA sample were performed for qRT-PCR analysis.We calculated the relative expression of each gene by the 2 −DDCt method (Livak & Schmittgen, 2001) and used the actins gene of S. caseolaris as a reference.

Physiological responses of S. apetala seedlings under chilling stress
After 3 h exposure to 4 • C, the leaves at the top of S. caseolaris seedlings began to show a slight curling, and after 6 h, the leaves in the top showed evident curling.After 12 h of chilling stress treatment, part of leaves had softened and exhibited loss of strength.After 24 h of 4 • C chilling stress, most of leaves had withered (Fig. 1).Electrolyte leakage, MDA content, soluble sugar content and soluble protein content increased significantly (Fig. 2).
The electrolyte leakage increased from 6.54 ± 0.46%(0 h) to 28.56 ± 5.32% (24 h); the MDA content increased from 2.15 ± 0.79 µg/g mol/g (0 h) to 14.01 ± 1.58 µg/g mol/g (24 h); the soluble sugar contents increased from 2.79 ± 0.76 µg/g mol/g (0 h) to 7.87 ± 1.67 µg/g mol/g(24 h); and the soluble protein content increased from 6.00 ± 0.84 µg/g mol/g (0 h) to 33.83 ± 2.27 µg/g mol/g (24 h).All of these parameters showed a linear increase.When exposed to chilling stress for 6 h, these parameters increased 118.68%, 299.98%, 160.67% and 56.15% compared with the control group.On the basis of these results, we selected the leaf samples exposed to 0 h and 6 h of chilling stress for RNA-seq to investigate the changes in genome-wide gene expression mechanism under chilling stress.The 0 h time the control (CK) and 6 h was regarded as the treatment group (CT), with three biological replicates.

Sequencing and de novo assembly of transcriptome
In total, we obtained 39.33 Gb clean bases from 6 cDNA libraries, and we obtained a total of 262.21 M clean reads from 280.46 M raw reads (Table S2).The clean reads were de novo assembled by Trinity software and they yielded 168,473 unigenes.All unigenes were longer than 200 bp and ranged from 297 to 16,287 bp.The 200-300 bp unigenes accounted for the largest proportion (13,159, 12.26%) (Fig. S1A).The mean unigene length was 1,762 bp.The N50 and N90 lengths were 2,849 bp and 937 bp respectively, and the GC content was 43.94%.

GO term enrichment analysis of DEGs
We identified a total of 3,706 DEGs, including 2,501 that were up-regulated and 1,205 that were down-regulated under the chilling stress (Fig. S2).The biological function information for the DEGs was analyzed by GO enrichment.A total of 2,902 DEGs were annotated in 1,411 terms and a total of 125 GO terms were identified.The most enriched GO categories among these DEGs were ''nucleus'', ''DNA binding'', ''metal ion binding'', ''DNA-binding transcription factor activity'', and ''regulation of transcription, DNA-templated '', with 366, 239, 204, 128, and 125 DEGs respectively (Fig. 3A).

Pathway enrichment analysis of DEGs
To characterize the active biological pathways of the DEGs in S. caseolaris under chilling stress, we classified the DEGs in KEGG biological pathways and conducted enrichment analysis.There were 2,557 DEGs annotated in 125 pathways.Among them, a total of nine KEGG pathways were significantly enriched with a P-value ≤ 0.05.The three most highly represented pathways were ''Plant hormone signal transduction (ko04075)'', ''Plantpathogen interaction (ko04626)'', and ''Glycerophospholipid metabolism (ko00564)'', with 115, 88 and 41 DEGs respectively (Fig. 3B).The largest number of DEGs included ''plant hormone signal transduction'', indicating that plant hormones play an important role in response to chilling stress in S. caseolaris.The second largest number of DEGs included ''Plant-pathogen interaction'', which suggested that low temperature exposed plants are more vulnerable to be pathogen infection.The ''Glycerophospholipid metabolism'' category had the third largest number of DEGs, suggesting that fatty acid metabolism responded positively to low temperature.A considerable number of biological reactions were altered in S. caseolaris in response to chilling stress.

Phytohormone related DEGs
KEGG analysis of DEGs in S. caseolaris showed that hormone signal transduction was significantly influenced by chilling stress (Fig. 3).There were 115 DEGs enriched in ''plant hormone signal transduction'' and most of them were up-regulated (Fig. 5).Especially in  In this GO analysis of DEGs, 41 DEGs enriched in the ''auxin-activated signaling pathway (GO:0009734)'' term.Among them, AUX1 (Unigene13181_All) is a key gene that encodes an auxin receptor, and down-regulated 2.50 fold.In the ABA signal pathway, PLY and PP2C are key enzymes controlled by the ABA signaling pathway.In this study, the genes encoding these two enzymes (CL1229.Contig3_All, Unigene14664_All, and CL5508.Contig1_All) were all up-regulated.Additionally, seven DEGs were enriched in the ''abscisic acid-activated signaling pathway (GO:0009738)'' term.In the ET signal pathway, we found that the one gene (CL2232.Contig4_All) encoding ETR down-regulated 6.81 fold.In the JA signal pathway, cetyl-CoA acyltransferase 1 is a key enzyme for synthesis of JA-CoA.In this study, the genes (CL10456.Contig6_All) encoding this enzyme were up-regulated by 3.93-fold.

ROS related DEGs
We found DEGs related to catalase, glutathione peroxidase, peroxidase, and superoxide dismutase were involved in chilling stress (Fig. 4B).Two genes encode superoxide

Validation of transcriptome data by RT-qPCR analyses
To validate the RNA-Seq results, we selected a total of 16 candidate genes for qRT-PCR analyses (Fig. 6).Based on enrichment analyses, all of these genes were involved in the chilling response.They represented different functional categories or pathways, such as sugar metabolism, antioxidant enzyme, phytohormones and TFs.The results indicated that the transcriptome results were reliable.

DISCUSSION
With its rapid growth adaptability, S. caseolaris is a good candidate species for mangrove afforestation in China.Low temperature, however, have limited its distribution in high latitude areas.Understanding the physiological and molecular mechanisms of chilling resistance can aid in breeding cold tolerant varieties of S. caseolaris.We performed physiological and transcriptomic analyses of leaves under 4 • C chilling stress, and then analyzed the DEGs involved in the antioxidant defense system, carbohydrate metabolism, plant hormone signals, and TFs.

Membrane injury of S. caseolaris seedlings caused by chilling stress
Membrane lipids are the first place that plants perceive and respond to abiotic stress signals such as low temperature.Lipids are involved in biological process, such as signal transduction, energy conversion, and metabolic regulation (Xie et al., 2018).Damaged membrane lipids increase electrolyte leakage, and a high electrolyte leakage rate reflects serious damage caused by chilling stress.Membrane permeability is an important physiological index of plant resistance to stress.MDA is the final degradation product of lipid peroxidation, and the MDA level reflects the damage of membrane lipid peroxidation and the stress tolerance of plant cells (Mutlu, Karadaolu & Nalbantolu, 2013;Wise & Naylor, 1987).Physiological analysis showed that the relative electrolyte conductivity and the MDA level both increased significantly after chilling stress.Membrane lipid peroxidation likely occurred and membrane lipids were damaged.Similar results were found in Magnolia wufengensis (Deng et al., 2019).Serious member lipid injury can result in seedling death.This might explain why most seedlings died when exposed to extreme low temperatures (Chen et al., 2017).We found that DEGs were enriched in ''Glycerophospholipid metabolism (ko00062)'' and ''Fatty acid elongation (ko00564)''.
The genes responded to the increased level of unsaturated fatty acids and membrane fluidity in the cells of S. caseolaris seedlings under chilling stress, which may be a mechanism of cellular self-repair.

Carbohydrate metabolism-related DEGs
When exposed to low temperatures, plants usually break down starch and increase the content of carbohydrates to maintain stability of cell membrane structure, this could increase the low temperature tolerance of plant cells (Guy et al., 2008).In GO enrichment analysis, 12 DEGs were enriched in ''starch binding (GO:2001070)'' (Fig. 3B).This could result in the conversion of starch into soluble sugar under chilling stress.Physiological experiments confirmed that the soluble sugars content of S. caseolaris increased after chilling stress.Increased levels of sucrose, raffinose, glucose, and fructose are involved in the low-temperature resistance of other plants (Anchordoguy et al., 1987;Ma et al., 2009).
The RNA-seq data showed that many DEGs that are involved in starch, fructose, sucrose, glucose, and trehalose were up-regulated in chilling stress treatment.These data were consistent with previous studies (Wang et al., 2018;Deng et al., 2019;Yang et al., 2019).
Plants can accumulate trehalose for oxidative detoxification to counteract abiotic stresses (Lai et al., 2017).In the DEGs, four genes were involved in the trehalose biosynthetic process and all of them were up-regulated, which is similar to the report of Deng et al. (2019).Sucrose synthase and sucrose phosphate synthase play important roles in sucrose metabolism and biosynthesis respectively (Winter & Huber, 2000).We found that two sucrose synthase genes and one sucrose-phosphate synthase gene were up-regulated.These genes may have been involved in the response to chilling stress.A similar result was found by Yang et al. (2019).Carbohydrate also regulate expression of jasmonate, ABA and COR genes (Rutten & Santarius, 1992;Welling & Palva, 2006).We found that many DEGs were involved in carbohydrate metabolism, which indicated that carbohydrates play an important role in S. caseolaris response to chilling stress.

Antioxidant defense system
Under environment stress, ROS such as superoxide, hydroxyl radicals and hydrogen peroxide can accumulate in plant cells and cause oxidative damage (Elsheery & Cao, 2008;Suzuki & Mittler, 2006).To scavenge ROS, antioxidant systems such as superoxide dismutase (SOD), catalase (CAT), peroxidase (POD), and ascorbate peroxidase (APX) are activated in response to plant stress (Mittler, 2002).SOD is the first antioxidant enzyme that could convert O 2− to hydrogen peroxide (H 2 O 2 ), and then POD and CAT could catalyze H 2 O 2 and convert it to H 2 O.We found DEGs related to CAT, glutathione peroxidase, POD, and SOD were up-regulated, indicating that they may have played a key role in activating antioxidant enzymes in S. caseolaris exposed to low temperatures.ROS are also signal molecules that can induce gene expression and protein synthesis to increase stress tolerance in plants (Jaspers & Kangasjärvi, 2010).ROS also is active the mitogen-activated protein kinase (MAPK) cascade pathway, which regulates low-temperature response genes involved in TFs and some hormones (Pitzschke, Schikora & Hirt, 2009;Xie et al., 2009).We found 15 DEGs enriched in ''MAP kinase activity'' and most of them were un-regulated, which indicated that MAPK signal transduction pathway also played an important role in chilling stress responses.These results indicated that many antioxidant genes are activated, and these would not only increase the activity of antioxidant enzymes to scavenge with ROS, but also active a sequence of signal to response chilling stress.

Phytohormone signal network responses to chilling stress
Phytohormones are the secondary signals in plant cells that play a significant role in adaptation to environmental stress (Peleg & Blumwald, 2011).They can activate a sequence of signal events and eventually induce the genes that directly respond to stress (Bari & Jones, 2009;Hakeem, Rehman & Tahir, 2014).We found four key response hormone signal transduction pathways (auxin, ABA, ET, and JA) that were significantly influenced by chilling stress.Auxin is a plant growth-promoting hormone (Nemhauser, Hong & Chory, 2006).In this study, 41 DEGs were enriched in the ''auxin-activated signaling pathway (GO:0009734)'' term.Down-regulation of the AUX1 gene may reduce plant growth in response to chilling stress.The ABA signaling pathway is important in regulating abiotic stress responses in plants (Danquah et al., 2014).It can regulate stomatal closure, reduce water loss and maintain cell growth (Peleg & Blumwald, 2011).PLY and PP2C are key genes controlled by the ABA signaling pathway.In present study, the genes were up-regulated in response to chilling stress.This could activate high expression of PP2C and promote stomatal closure and decrease transpiration to reduce water loss.Similar results have been found on Arabidopsis thaliana (Tähtiharju & Palva, 2001) and Oryza sativa L (Guan et al., 2019).Additionally, the DEGs were also enriched in the ''abscisic acid-activated signaling pathway (GO:0009738)'' term, which further confirmed that ABA-dependent signaling is involved in the low temperature response of S. caseolaris.Genes related to ET biosynthesis often down-regulated expressed under low temperature treatment (Chu & Lee, 1989).In this study, the key gene encoding ETR was down-regulated, which was consistent with previous reports (Yang & Huang, 2018).JA plays a significant role in plants responding to abiotic and biotic stress (Wasternack & Song, 2017;Wasternack & Strnad, 2018).The up-regulated key genes may promote the synthesis of JA-CoA and promote the synthesis of JA.Exogenous methyl salicylic acid (MeSA) or methyl jasmonate (MeJA) could induce the transcription of pathogenesis-related genes, and increase chilling tolerance and pathogen resistance (Ding et al., 2002).This may explain why many DEGs were enriched in ''Plantpathogen interaction''.Similar results have been found in K. obovata, M. wufengensis, and Camellia sinensis under chilling stress (Deng et al., 2019;Li et al., 2019;Su et al., 2019).However, how those plant hormones mediated signaling participates in chilling stress responses in S. caseolaris still need to be further study.

TFs responding to chilling stress
TFs play important roles in plant growth and stress tolerance (Lee, Henderson & Zhu, 2005).When a plant experiences the chilling stress, a series of signal transduction pathways could activate TFs.The activated TFs will activate the expressions of a cascade of downstream resistance-related genes to enhance the cold tolerance by special binding (Kasuga et al., 1999).Many TF families such as DREB (Ahmed et al., 2017;Chen et al., 2012), MYB (Zhao et al., 2012;Agarwal et al., 2006;Chen et al., 2006), andWRKY (Zou, Jiang &Yu, 2010) are related to plants resistance to low temperature stress.We identified 232 TFs that were classified into 41 families.AP2-ERF, MYB, C3H, WRKY, and ARF were five major TF families identified as DEGs, suggest that they played a crucial role in response to chilling stress.AP2-ERF genes belong to a plant-specific TF family involved in regulating stress responses (Mizoi, Shinozaki & Yamaguchi-Shinozaki, 2012).The AP2-ERF superfamily consists of the ERF, RAV, AP2, and DREB families.Overexpression of AP2 TF, JcDREB, and BpDREB2 increased cold tolerance of Arabidopsis (Tang et al., 2011;Sun et al., 2014).Previous studies verified that AP2-ERFs involved in a lot of plant hormone mediated abiotic stress responses, such as ABA, ET, GA, CTK, and BR (Colebrook et al., 2014;Kazan, 2015;Tao et al., 2015;Nolan et al., 2017;Sah, Reddy & Li, 2016).MYB are involved in chilling stress and regulation of ABA-responsive genes, and constitute with ABA-dependent, induced and mediated types.Which plays an important role in the upstream stage of low temperature transduction (Gyoungju et al., 2016).AtMYB52, AtMYB70, and AtMYB82 were ABA-dependent genes in responding abiotic stresses of Arabidopsis (Park, Kang & Kim, 2011); TaMYB1, TaMYB2, TaMYB3R1 were ABA-induced genes in responding abiotic stresses of wheat (Cai et al., 2011;Tong et al., 2006); AtMYB2, AtMYB41 were ABA-mediated types that promote the accumulation of ABA in responding abiotic stresses (Abe et al., 2003;Lippold et al., 2009).The WRKY genes can respond strongly and rapidly to abiotic stresses (Chen et al., 2012).The WRKY46 gene enhanced the low temperature tolerance of transgenic cucumber by activating a series of cold-stress responsive genes in an ABA-dependent manner (Zhang et al., 2016).MaWRKY25 can be induced by low temperature or MeJA, which could increase the low temperature tolerance of banana.It could also active the MaLOX2, MaAOS3, MaOPR3 genes to promote the synthesis of JA by special binding.Which could reduce the damage of banana fruit by low temperature (Ye et al., 2016).Other TFs, such as the HSF, bZIP, and GRAS families, are induced by chilling stress, and members of these families can affect the cold resistance of plants (Chinnusamy, Zhu & Sunkar, 2010).These results indicated that many TFs were involved in the seedlings of S. caseolaris in response to chilling stress through different pathways.Besides, the both un-regulated TFs and phytohormones maybe have some connection, and the regulatory network need further study.

CONCLUSIONS
We analyzed at the morphological, physiological and transcriptome level of S. caseolaris response to chilling stress.Chilling stress of S. caseolaris likely involves damage to the

Figure 4
Figure 4 Heat map of relative expression levels of DEGs involved in main transcription factors, sugar metabolism and the antioxidant defense system.(A) The main sugar metabolism-related genes.(B) The main antioxidant-related genes.(C) The main transcription factors.Full-size DOI: 10.7717/peerj.11506/fig-4

Figure 7
Figure 7 Hypothetical model of the events occurring in the S. caseolaris leaves under chilling stress.Full-size DOI: 10.7717/peerj.11506/fig-7