Transcriptomic analysis identifies novel genes and pathways for salt stress responses in Suaeda salsa leaves

Salinity is a critical abiotic stress, which significantly impacts the agricultural yield worldwide. Identification of the molecular mechanisms underlying the salt tolerance in euhalophyte Suaeda salsa is conducive to the development of salt-resistant crops. In the present study, high-throughput RNA sequencing was performed after S. salsa leaves were exposed to 300 mM NaCl for 7 days, and 7,753 unigenes were identified as differently expressed genes (DEGs) in S. salsa, including 3,638 increased and 4,115 decreased unigenes. Moreover, hundreds of pathways were predicted to participate in salt stress response in S. salsa by Gene Ontology (GO), MapMan and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses, including ion transport and sequestration as well as photoprotection of photosystem (PS) II. The GO enrichment analysis indicated that genes related to ion transport, reactive oxygen species (ROS) scavenging and transcriptional factors were highly expressed upon NaCl treatment. The excessive Na+ and Cl− ions were supposed to be absorbed into the vacuole for ion sequestration and balance adjustment by potassium transporters (such as KEA3) with high expressions. Moreover, we predicted that mutiple candidate genes associated with photosynthesis (such as PSB33 and ABA4), ROS (such as TAU9 and PHI8) and transcriptional regulation (HB-7 and MYB78) pathways could mitigate salt stress-caused damage in S. salsa.

functional annotations of unigenes in S. salsa. Similarity test was carried out to annotate unigenes against different databases by BLASTX. All 86,255 (100%) unigenes were annotated in at least one database. A total of 64,872 (75.20%) (Online Resource 1), 49,254 (57.10%) and 48,351 (56.05%) unigenes showed demonstrated similarity to sequences in NT, NR and PFAM databases, respectively, with an E-value threshold of 1e −5 (Fig. 1A). Moreover, we assigned 48,463 (56.18%) unigenes, with an E-value cutoff of 1e −6 , in Gene Ontology (GO) database using Blast2GO v2.5. Besides, 59,738 unigenes of the S. salsa were assigned to A. thaliana gene IDs for GO annotation mapping by BLASTX with an E-value cutoff of 1e −5 and categorized into GO terms for GO analysis.  Table S1).

GO, MapMan and KEGG enrichment results of DEGs in S. salsa.
To uncover the molecular mechanisms underlying the salt tolerance of S. salsa leaves, the DEGs were characterized with GO databases. Consequently, 171 biological process (BP) terms were enriched by 3,638 up-regulated unigenes with the cutoff of P value < 0.05, including "calcium ion transmembrane transport" (GO: 0070588), "phosphate ion homeostasis" (GO: 0055062) and "cellular protein catabolic process" (GO: 0044257) (Tables 2; S2). The 4,115 genes with decreased expressions were identified to be enriched in 208 BP terms, such as "translation" (GO: 0006412), "polyamine catabolic process" (GO: 0006598) and "photosynthesis, light harvesting in photosystem I" (GO: 0009768) (Tables 2; S2). The DEGs were dispatched to 1,595 and 2,328 homologs in Arabidopsis, respectively. We assigned www.nature.com/scientificreports www.nature.com/scientificreports/ these genes to 905 pathways by MapMan. Among these pathways, the dysregulated genes were enriched in 55 pathways, with a cutoff of P value < 0.05.
Due to the large numbers and complex branch structures of GO categories, REVIGO was used to find representative subgroups of the terms using a simple clustering algorithm that relies on semantic similarity measures. The 171 BP terms enriched by genes with increased expressions were integrated into 10 groups ( Fig. 2A), 37 terms were dispatched to "ribonucleoside monophosphate biosynthesis" subset, and 32 terms were classified to "response to absence of light" group. Moreover, 208 BP terms enriched by genes with decreased expressions were dispatched into 10 subsets (Fig. 2B), such as "translation" (including 52 terms) and "cellular response to sulfate starvation" (including 27 terms). Figure 3 illustrates the metabolism result of MapMan analysis. The KEGG pathway "plant hormone signal transduction" (ko04075) was enriched by 23 up-regulated unigenes, of which the cell enlargement-related sub-pathway "tryptophan metabolism" was enriched by auxin-responsive protein (AUXIAA, Cluster-12522.63856, L 2 fc = 1.341), auxin response factor (ARF, Cluster-12522.18986, L 2 fc = 5.928) and SAUR family protein (SAUR, Cluster-1,512.0, L 2 fc = 6.226); the stomatal closure regulation sub-pathway "carotenoid biosynthesis" was enriched by four up-regulated unigenes, including abscisic acid receptor PYR/PYL family protein (PYR, Cluster-12522.12070, L 2 fc = 4.312), protein phosphatase (PP2C, Cluster-12522.20854, L 2 fc = 9.005) and ABA responsive element binding factor (ABF, Cluster-12522.47885, L 2 fc = 2.532).

The differentially expressed transcription factors (TFs).
In the present study, we found that 2,656 unigenes were dispatched to 47 TF families (Table S3), and 293 differentially expressed TFs were found in NaCl-treated leaf samples, including 145 genes with increased expressions and 148 gens with decreased expressions (Table S3). Of the up-regulated TFs, the largest number of up-regulated genes was found in HB (18 unigenes), followed by MYB (13 unigenes) and bZIP families (10 unigenes). Conversely, the largest number of down-regulated genes was found in MYB (20 unigenes), followed by the bHLH and AP2/ERF families, which contained 13 and 10 unigenes with decreased expressions, respectively (Table 3).

Real-time quantitative PCR (RT-qPCR) validation.
To verify the RNA-Seq results, an alternative strategy was selected for the dysregulated unigenes. Four up-regulated and four down-regulated unigenes were selected for RT-qPCR validation with the same RNA samples used for RNA-Seq. Primers spanning exon-exon junctions were designed (Table S4). Figure 4 reveals that these two methods demonstrated similar gene expression profiling in most cases. For example, the ortholog of salt stress responding gene ITN1, Cluster-12522.40169, was significantly increased by RT-qPCR method (Fig. 4), and it was also an up-regulated unigene in the salt-treated samples by RNA-Seq (L 2 fc = 1.573).

Discussion
Among the most serious environmental factors, salt stress greatly affects the development and yield of crops worldwide. External salt concentrations severely inhibit photosynthesis in glycophytes 1,5,57,58 , while halophytes maintain or have enhanced photosynthesis and ion balance at these salt concentrations 36,37,59 . S. salsa, a euhalophyte plant with salt-diluting ability, can adapt to high salinity, and its photosystem (PS) II shows high resistance to salt stress 41,44,[60][61][62][63] . Many previous studies have shown the responses of S. salsa under salt stress. However, the underlying mechanisms of salt deposition via ion compartmentalization, osmotic adjustment and transcriptional regulation remain largely unclear. Therefore, we attempted to explore the essential salt responding genes and their modulators in S. salsa. Our transcriptome analysis provided valuable insights into the molecular mechanisms of salt tolerance and insulation of photosynthesis in S. salsa leaves.
Salt stress induces damage primarily through osmotic and ion stresses. Therefore, halophytes can resist to salt stress due to their ability to adjust osmotic state and ion balance 14,28,64 . The excess Na + , Cl − and oxidative stress in the intracellular or extracellular environment activate the osmotic adjustment or homeostasis regulating salt stress responses 27,28 . In salt-diluting halophyte, the vacuole is supposed to absorb excessive Na + and Cl − for ion sequestration and balance adjustment. Na + generally enters cells via ion transporters, such as potassium transporter (KT) [27][28][29] . Guo et al. 65 have shown that the expression levels of Na + /H + exchanger, V-H + ATPase, choline monooxygenase, potassium and chloride channels are up-regulated in S. salsa leaves upon 500-mM saline treatment, which can reduce the over-accumulation of Na + and Cl − . Our transcriptomic data showed that the expressions of KT-encoding genes were increased to enhance protein yield to absorb more ions, such as K + efflux antiporter 3 (Cluster-12522.45367, L 2 fc = 7.991). Furthermore, calcium ion transmembrane transport participating cation exchanger 3 (Cluster-12522.69849, L 2 fc = 6.425) and 5 (Cluster-12522.41534, L 2 fc = 3.393) were found to be highly expressed under salt stress, indicating their functions in transporting ions in S. salsa. Besides, excess Cl-was predicted to be transported by chloride channel B (Cluster-12522.3333, L 2 fc = 3.546) and C (Cluster-12522.30526, L 2 fc = 3.828).
In plants, salt stress-triggered ion injury and osmotic stress usually affect a set of basic biosynthetic functions 27,28 . As one of the most sensitive physiological indices, photosynthesis is susceptible to salt stress 7,44,49,50 . Although S. salsa shows significant resistance not only to salinity stress but also to photoinhibition even when exposed to 400 mM NaCl or full sunlight 41,60,61 , the photosynthesis rate is decreased by Cl − treatments 66     The representatives are joined into "superclusters" of loosely related terms, visualized with different colors. Size of the rectangles was adjusted to reflect the P value of the GO term calculated by TopGO. In this study, 171 upregulated processes were integrated into 10 groups, 13 subsets were summarized for 208 biological processes enriched by down-regulated genes. (2020) 10:4236 | https://doi.org/10.1038/s41598-020-61204-x www.nature.com/scientificreports www.nature.com/scientificreports/ PSB33 (Cluster-12522.32409, L 2 fc = 9.833) was up-regulated in NaCl-treated leaves in S. salsa. Besides, neoxanthin acts as an antioxidant within the PSII supercomplex in thylakoids, involving in the photoprotection of PSII 67 . In this study, the expressions of neoxanthin biosynthesis-required gene abscisic acid-deficient 4 (ABA4, Figure 3. Global view of differently expressed genes (DEGs) involved in diverse metabolic pathways. DEGs genes were selected for the metabolic pathway analysis using the MapMan software (3.5.1 R2). The colored boxes indicate the Log 2 of expression ratio of DEGs genes. The dys-regulated unigenes were assigned to 1,595 and 2,328 homologs in Arabidopsis, respectively. These genes were mapped to 905 pathways by MapMan, of which, 55 pathways were filtered enriched by the dys-regulated genes with the cutoff p-value <0.05.

TF family
Up gene num. Down gene num.  www.nature.com/scientificreports www.nature.com/scientificreports/ Cluster-12522.63197, L 2 fc = 1.103) was increased when exposed to salt stress. Our transcriptomic data provided solid evidence supporting the tolerance of PSII to salt stress for euhalophyte S. salsa 68 .

MYB
Saline stress inevitably results in the ROS overproduction in plants. Therefore, the activation of the ROS scavenging system is necessary for the plant to avoid additional oxidative stress 69,70 . An effective enzymatic antioxidant defense system has been evolved by plants for scavenging ROS [71][72][73] . In this study, the expressions of ascorbate peroxidase (APX), glutathione S-transferase (GST) and glutathione peroxidase (GPX) were increased to varying degrees to avoid oxidative stress induced by salinity, including ascorbate peroxidase 1 (Cluster-12522.1335, L 2 fc = 9.348), glutathione S-transferase TAU9 (Cluster-12522.44190, L 2 fc = 4.923), glutathione S-transferase PHI8 (Cluster-12522.25968, L 2 fc = 8.422) and glutathione peroxidase 4 (Cluster-12522.20682, L 2 fc = 8.949). Our transcriptomic data showed that S. salsa leaves could resist the oxidative stress caused by high salinity by recruiting the ROS scavenging system.
TFs play critical roles in salt tolerance in plants via transcriptional regulation of the downstream genes responsible for response to salt challenges. Members of several TF families were verified to control and modulate salt stress adaptive pathways, as bZIP, WRKY, ERF/AP2, HB, MYB and bHLH [74][75][76][77] . The HB family TFs with homeodomain-leucine zipper domain are transcriptionally regulated in an ABA-dependent manner and may act in a signal transduction pathway which mediates abiotic stress response. In the present study, 18 and one HB TFs were up-regulated and down-regulated in NaCl-treated samples, respectively, including homeobox 7 (Cluster-12522.23710, L 2 fc = 2.558) and BEL1-like homeodomain 2 (Cluster-12522.57518, L 2 fc = 3.612). The MYB family members have been proved to be key factors in regulatory networks of abiotic stresses 75,77 . In our current study, 13 and 20 MYB family TFs were up-regulated and down-regulated in NaCl-treated samples, respectively, including MYB domain protein 78 (Cluster-10,720.0, L 2 fc = 3.050) and Reveille 8 (Cluster-12522.18843, L 2 fc = 3.777). Numbers of bZIP proteins isolated from diverse species, including Arabidopsis, rice, tomato and halophytes, can enhance the salt tolerance ability 77 . Moreover, 10 bZIP family members were up-regulated in salt-treated samples in this analysis, such as unfertilized embryo sac 12 (Cluster-12522.12729, L 2 fc = 8.992) and iaa-leucine resistant 3 (Cluster-12,694.1, L 2 fc = 2.701). The spatio-temporal gene expression patterns of these TFs and novel genes involved in RNA processing in S. salsa could be helpful in understanding the transcriptional regulation in euhalophytes under salt stress.
Plant hormones, especially auxin, ethylene and jasmonic acid signaling transduction pathways, were all up-regulated in response to saline treatment in S. salsa, which are important to gene regulations of ion transport and antioxidation 65 . It has been well documented that as an endogenous signaling molecule, ABA enables plants to survive under salt stress 78,79 . In the present study, we investigated the ABA signaling pathway in response to saline stress in the leaves of S. salsa, and found that a part of molecular mechanism responding to NaCl could be attributed to the ABA-dependent signaling pathway. The osmotic stress-induced ABA signaling component aldehyde dehydrogenase 7B4 (Cluster-12522.47012, L 2 fc = 7.613) and aldehyde dehydrogenase 3H1 (Cluster-12522.7455, L 2 fc = 1.033) 80 are dramatically activated in NaCl-treated leaves. Besides, the expressions of MYB/NAC TF family members were increased when exposed to salt stress in this study, including ABA-induced MYB domain protein 78 (Cluster-10,720.0, L 2 fc = 3.050) and NAC domain-containing protein 83 (Cluster-12522.39867, L 2 fc = 1.533). Thus, ABA might exert direct effects on MYB/NAC TFs and modulate salt-resistance in S. salsa under salt stress.

conclusions
In the present study, a transcriptome analysis was carried out to assess the salt stress responses in S. salsa leaves. RNA-Seq revealed that the expressions of 3,638 unigenes were increased when exposed to salt stress, while 4,115 unigenes were down-regulated in S. salsa. Moreover, GO, MapMan and KEGG enrichment analyses supported that hundreds of pathways were involved in salt stress response in S. salsa, such as ion transport and sequestration as well as photoprotection of PSII. The GO enrichment analysis indicated that genes associated with ion transport, ROS scavenging and TFs were highly expressed under salt stress. The excessive Na + and Cl − were transported into the vacuole for ion sequestration and balance adjustment by highly expressed potassium transporters (such as KEA3). In addition, various candidate genes associated with photosynthesis (such as PSB33 and ABA4), ROS (such as TAU9 and PHI8) and transcriptional regulation (HB-7 and MYB78) pathways were predicted to ameliorate salt stress-induced damage in S. salsa.

Materials and Methods
Plant material and stress treatments. Seeds of S. salsa were collected at the Yellow River Delta, Dongying, China (37 °7′38″N; 118 °41′28″E), in October 2017. The plump seeds were firstly screened, disinfected with 75% alcohol for 2 min, and then washed by water. After vernalization at low temperature (4 °C) for 3 days and germination for 2 days, the seedlings of S. salsa were grown in plastic pots (14 cm in diameter and 13 cm in height) filled with sands in a greenhouse under the conditions of 15/9-h light/dark cycle at temperature (28 ± 3/20 ± 3 °C, day/night) with illumination of 600 μmol m −2 s −1 , and half-strength Hoagland nutrient solution was given daily.
No physiological interpretation of the data was included in this study because the design of this study referenced several existing studies demonstrating the changes/regulations of S. salsa under salt stress at the morphological and physiological levels. Lu et al. 60 and Sui et al. 81 have demonstrated the physiological meaning of salt treatment/absence for S. salsa. Compared with the control (0 mM NaCl treatment) group, the growth of S. salsa is significantly increased by salt treatment. Sui et al. 81 have shown that seedlings with 2-3 branches are subjected to 300 mM NaCl treatment, and the growth of S. salsa is significantly increased. PSII is rather tolerant to NaCl, and salt stress has no effect on PSII photochemistry in dark-adapted leaves. Photosynthetic electron transport is increased under salt stress. Although 400-500 mM NaCl is regarded as high salinity for S. salsa, 300 mM NaCl is the inflection point of the growth, ion accumulation and enzymatic antioxidant defense system for scavenging of ROS 81 . Given these above-mentioned physiological evidence, 300 mM NaCl was selected for salt stress treatment.
Vigorous seedlings were selected and transferred into plastic pots containing well-washed river sand (three plants per pot) when the sixth leaf emerged. NaCl was dissolved in the full-strength Hoagland nutrient solution Plants of control group was irrigated with full-strength Hoagland nutrient solution (0 mM NaCl), while others were exposed to 300 mM NaCl. The NaCl level was gradually increased by 50 mM per day to avoid osmotic shock. After 7 days, the seedlings were selected to determine the indices. Three replicates of each treatment were sampled for transcriptome sequencing.
Illumina library construction and sequencing. Sequencing libraries were generated using the NEBNext ® Ultra TM RNA Library Prep Kit for Illumina ® (NEB, USA) based on the manufacturer's instructions, and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. The first-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H − ). Second-strand cDNA was synthesized using DNA polymerase I and RNase H. The AMPure XP system (Beckman Coulter, Beverly, USA) was employed to purify cDNA fragments of 150~200-bp. The size-selected, adaptor-ligated fragments were purified and enriched by PCR amplification. The resulting products were used for sequencing analysis. The Illumina HiSeq X platform (Illumina, San Diego, CA, USA) was used to perform high-throughput sequencing. The full assembly data were submitted to the NCBI Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra), SRA accession: PRJNA527358. The processed data files including the assembled sequences and abundance measurements were uploaded to the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), GEO accession: GSE145366.
De novo assembly of the S. salsa transcriptome. RNA sequencing and de novo transcriptome assembly were conducted to create reference sequence libraries for S. salsa. For each species, an RNA sample of every accession was separately sequenced. cDNA library construction and Illumina paired-end 150 bp sequencing (PE150) were performed at Novogene Co., LTD. (http://www.novogene.com/) according to instructions provided by Illumina Inc. Raw reads of FASTQ format were firstly processed through in-house perl scripts. In this step, clean reads were obtained by removing adapter-containing reads, ploy-N-containing reads and low-quality reads from raw data. Meanwhile, Q20, Q30, GC content and sequence duplication level of the clean data were calculated. All the downstream analyses were carried out based on high-quality clean data. The remaining high-quality reads were used for transcriptome assembly using the Trinity software pipeline with default parameters 82 . Clean data sets of three accessions were pooled for de novo assembly and comprehensive sequence library construction in each studied species. De novo assembled unigene sequences were used for BLAST searches and annotation against public databases (NR, NT, Swiss-Prot, Pfam, KOG/COG, Swiss-Prot, KEGG Ortholog database and Gene Ontology) with an E-value threshold of 1e −5 .

Scientific RepoRtS |
(2020) 10:4236 | https://doi.org/10.1038/s41598-020-61204-x www.nature.com/scientificreports www.nature.com/scientificreports/ Determination of gene expression in S. salsa. Six independent transcript libraries were generated for three accessions of S. salsa by a PE150 sequencing analysis. Gene expression levels for all samples were estimated by RSEM 83 . The clean reads were aligned to the de novo assembled transcriptome. Gene expression levels in leaf samples of S. salsa were calculated by the FPKM method 84 . The expression levels between salt-treated and control samples were compared according to the FPKM values, with a cutoff of Padj <0.05 and L 2 fc >1.

GO, MapMan and KEGG annotation and enrichment.
The unigenes of S. salsa were mapped to A. thaliana gene IDs by sequence similarity searching against the genome of A. thaliana with an E-value cutoff of 1e −5 . The topGO package of R 85 was used to perform the GO enrichment analysis for DEGs in NaCl-treated samples. The DEGs in S. salsa were annotated onto metabolic pathways using MapMan (version 3.5.1 R2) 86 . The DEGs of S. salsa unigene IDs were transferred to the Arabidopsis TAIR locus IDs during the MapMan analysis. The software KOBAS was used to test the statistical enrichment of DEGs in KEGG pathways in S. salsa 87 .
RT-qPCR verification. The expression patterns revealed by the RNA-Seq analysis were validated using the RT-qPCR. The extracted RNA samples were treated with DNaseI and reversely transcribed into cDNA using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China). Five up-regulated and five down-regulated unigenes in S. salsa were selected for the RT-qPCR assay, including Cluster-12522.  (Table S4). qPCR was performed using SYBR Green qPCR Master Mix (DBI, Germany) on an ABI7500 Real-Time PCR System (ABI, USA). Three replicates were performed for each experiment, and the amplification specificity was evaluated by melting curve analysis. Relative expressions of target genes were calculated using the 2 −ΔΔCt method. All values were presented as the mean ± standard deviation (SD). Significant differences were determined using GraphPad software (version: 5.0). The t-test was used to analyze the relationship between the expressions of control and NaCl-treated samples.

Consent for publication.
We have carefully read and adhered to editorial policies for manuscripts. We declare that the content of this manuscript has not been published or submitted for publication elsewhere.

Data availability
We declare that the all data and materials of this manuscript including the supplementary datasets are available to the journal and all readers.