De novo transcriptome sequencing and analysis of genes related to salt stress response in Glehnia littoralis

Introduction

High soil salinity has caused extremely negative effects on global agricultural production and ecological environments. Plant growth, development, metabolism, and protein functions can be retarded and inhibited under salt stress (Carillo et al., 2011). High salinity affects plants during two phases, the early and rapid osmotic stress phase and the later and slower ionic toxicity phase (Munns & Tester, 2008). To resist damage caused by salt stress and survive in high-salt conditions, plants have gradually formed a series of complex mechanisms including osmotic regulation, Na⁺ exclusion, and Na⁺ compartmentalisation (Deinlein et al., 2014; Munns & Tester, 2008). A variety of signal transduction pathways are involved in plant responses to salt stress, including the phospholipid signaling pathway (Hong et al., 2016), salt overly sensitive (SOS) pathway (Qiu et al., 2002; Quan et al., 2007), abscisic acid (ABA) pathway (Hauser et al., 2017), calcium-dependent protein kinase (CDPK) pathway, and mitogen-activated protein kinase (MAPK) cascade pathway (Zhu, 2016). These pathways interact to form signal transduction networks in plant responses to salt stress. Numerous breakthroughs have been reported in the study of plant salt tolerance; however, many salt-related functional genes remain to be discovered and identified.

Glehnia littoralis Fr. Schmidt ex Miq. belongs to the Umbelliferae family. Glehnia littoralis germplasm resources are distributed in Japan, Russia, China, the Korean Peninsula, the Pacific Rim, and coastal areas of the United States. It is an endangered halophyte shown to exhibit salt tolerance (Peng et al., 2014), and has traditional and medicinal uses. Peeled and dried G. littoralis roots, named Radix Glehniae, have been used as a traditional Chinese herbal medicine for moistening lungs, removing phlegm, curing respiratory and gastrointestinal disorders, and anti-inflammation. This herbal medicine is also a major tonic component in anti-aging and health promotion prescriptions. Coumarins and polyacetylenes are the primary active constituents of G. littoralis (Yuan et al., 2002). Although G. littoralis has great ecological and commercial value, salt-related mechanisms in G. littoralis remain largely unknown. Thus far, most studies of salt tolerance in G. littoralis have focused on anatomical and morphological adaptations to high-salinity environments (Voronkova et al., 2011); for example, Glehnia littoralis leaves exhibit dorsoventral structure and have developed secretory trichomes that remove excess salts. However, the salt tolerance mechanisms and functional gene information of G. littoralis have not been explored.

In recent years, de novo transcriptome assembly and annotation of non-model plant species, especially those with unknown genomic sequences, have become possible through RNA sequencing (RNA-Seq) technology (Ward, Ponnala & Weber, 2012; Xiao et al., 2013). Transcriptome analysis has greatly facilitated the discovery of putative functional genes or proteins involved in various plant biological processes (Cao & Deng, 2017; Lu et al., 2017; Wang et al., 2017). Comprehensive transcriptome analysis has been performed on various plant species in response to various abiotic stresses, such as cold, drought, salt, and hormone stresses, and has provided insight into their functional gene expression and regulation, signal transduction networks, and metabolite biosynthetic pathways (Postnikova, Shao & Nemchinov, 2013; Wang et al., 2013; Wang et al., 2017; Zhang et al., 2014a; Zhou et al., 2017). For example, recently in Litchi chinensis, 73,117 unigenes were assembled and 11,741 unigenes were identified as being both chilling and ROS responsive genes (CRRGs) (Lu et al., 2017). Some of these CRRGs are involved in flowering, plant hormone signal transduction and plant hormone biosynthesis, and exhibit relationships within genes co-expression networks. Zhang et al. (2014a) performed transcriptome analysis of Populus pruinosa in response to salt stress by constructing six libraries from control and salt-treatment calli at different time-points. They identified 9,216 differentially expressed genes (DEGs) and found that during the salt treatment, most of the DEGs were activated early (within 24 h) and stabilized after 48 h.

Halophytes are ideal plants for studying plant salt tolerance genetics using transcriptome sequencing. Comparative transcriptome analyses have previously been reported for a few halophytes, such as Halogeton glomeratus (Yao et al., 2018; Wang et al., 2015), Iris lactea var. chinensis (Gu et al., 2018), Kochia sieversiana (Zhao et al., 2017), Nitraria sibirica Pall. (Li et al., 2017a), Suaeda maritima (Gharat et al., 2016), Sporobolus virginicus (Yamamoto et al., 2015), and Mesembryanthemum crystallinum (Oh et al., 2015). In one such study, Wang et al. (2015) compared the transcriptomes of H. glomeratus exposed to NaCl for 6, 12, 24, and 72 h and discovered 2,223, 5,643, 7,510 and 10,908 DEGs, respectively. Later, Yao et al. (2018) analysed root transcriptomes of H. glomeratus under salt stress and identified the core DEGs regulating Na⁺ uptake and transport in the roots. G. littoralis is a typical halophyte belonging to the family Umbelliferae. To date, no genomic or transcriptome sequencing information about salt stress responses has been reported for any halophytes belonging to Umbelliferara, including G. littoralis. As a medicinal halophyte, salt-associated molecular mechanisms in G. littoralis are important to a comprehensive understanding of this plant. In the current study, we performed comparative transcriptome sequencing of G. littoralis under NaCl treatment, and identified a total of 105,875 unigenes, 71.37% of which were annotated in public databases. These transcriptome analysis results will help to promote genomics and genetics studies on salt-associated mechanisms and on secondary metabolite biosynthesis in G. littoralis.

Materials & Methods

Results

G. littoralis transcriptome sequencing profile under NaCl treatment

To perform a comprehensive analysis of the G. littoralis transcriptome under salt treatment, we constructed two RNA-Seq libraries, representing non-salt treatment (0 h, GL1) and 24 h of 200 mM NaCl treatment (GL2), using the Illumina HiSeq 2000 platform. The GL1 and GL2 libraries produced 5.31 and 5.09 Gb of clean data from the paired-end reads with Q20 percentages of 96.7% and 96.5%, respectively (Table 1). The clean reads from the GL1 and GL2 libraries were used for de novo assembly with the Trinity package. The Trinity method determines splice isoforms, distinguishes transcripts from recent duplicates, and identifies allelic variants (Grabherr et al., 2011); thus, each isoform corresponded to a unigene in the present dataset. In total, there were 105,875 unigenes (isoforms) of between 351 and 15,667 bp in sequence assembly, with a mean length of 1,314.44 bp (Fig. 1, Table 2).

Table 1:

Summary of G. littoralis RNA sequencing (RNA-Seq).

Samples	Raw reads	Clean reads	Clean bases	Q20 percentage (%)
GL1 (control)	54568070	53715849	5.31 Gbyte	96.70%
GL2 (NaCl)	52469596	51566060	5.09 Gbyte	96.50%

DOI: 10.7717/peerj.5681/table-1

Table 2:

Summary of G. littoralis transcriptome assembly.

Type	Number
Total genes	53,092
Total unigenes	105,875
Total residues	1.39E+08
Average length	1,314.44
Largest unigene	1,5667
Smallest unigene	351

DOI: 10.7717/peerj.5681/table-2

Sequence annotation and classification

To annotate the total unigenes, sequences were aligned using public databases (Table 3). In total, 75,535 unigenes (71.34%) were matched in the NR database, 25,871 unigenes (24.44%) in the KOG database, 47,043 unigenes (44.43%) in the GO database, and 24,864 (23.48%) in the KEGG database (Table 3). Taken together, 75,559 unigenes (71.37%) were successfully annotated in at least one of the databases (Table 3). In the GO annotation analysis, 47,043 unigenes were found and divided into three groups. The proportions of unigenes aligned in each GO term are shown in Fig. 2. Based on sequence homology with genes with known functions, these unigenes could be assigned to one or more ontologies, including 132,455 unigenes assigned to the biological process category, 100,205 unigenes assigned to the cellular component category, and 58,269 unigenes assigned to the molecular function category (Table S2). Some unigenes were assigned to more than one GO term, that is, a unigene may belong to several GO term notes and each GO term may correspond to multiple genes. Cellular process (28,273) and metabolic process (27,405) were highly represented groups within the biological process category. In the cellular component category, cell (23,485), cell part (23,483), and organelle (17,311) were highly represented groups. Moreover, catalytic activity (26,412) and binding (24,277) constituted the largest proportion in the molecular function category (Fig. 2, Table S2).

Table 3:

Functional annotations of all unigenes found in the public databases.

Annotation	Number of unigenes	Percentage (%)
Annotated in NR	75,535	71.34356553
Annotated in KOG	25,871	24.43541913
Annotated in GO	47,043	44.4325856
Annotated in KEGG	24,864	23.48429752
Annotated in at least one database	75,559	71.36623377
Total unigenes	105,875	100

DOI: 10.7717/peerj.5681/table-3

In addition, 25,871 unigenes had significant matches in the KOG database for functional prediction and classification (Fig. 3, Table 2). Among the 25 KOG categories, the cluster for general function prediction only (5,358) was the largest category, followed by signal transduction mechanisms (3,494), posttranslational modification, protein turnover, chaperones (2,630), and transcription (1,766). Cell motility (five) and extracellular structures (70) were the smallest categories (Fig. 3, Table S3). We also performed KEGG pathway analysis of the unigenes. In total, 41,556 unigenes participated in 25 predicted metabolic pathways (Fig. 4; Table S4), which were divided into five categories, including metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems. The largest category was metabolism, including global and overview maps (10,892), followed by carbohydrate metabolism (3,727), amino acid metabolism, (1,987), lipid metabolism (1,931), energy metabolism (1,500), nucleotide metabolism (902), glycan biosynthesis and metabolism (703), metabolism of terpenoids and polyketides (678), metabolism of other amino acids (670), metabolism of cofactors and vitamins (663), biosynthesis of other secondary metabolites (575), and xenobiotics biodegradation and metabolism (539) (Fig. 4, Table S4).

Identification of DEGs under NaCl treatment

The clean reads from GL1 (control) and GL2 (NaCl) were mapped to the assembled transcriptome sequence to acquire read counts data using RSEM. Then, DEGs between two samples were identified using edgeR package with FDR <0.05 and —log₂ fold-change— ≥ 1. Unigene expression was calculated using the mappable reads, and the results were normalised to FPKM values. For each sample, more than 99% of FPKM values were between 1 and 100 (Fig. 5A). In total, we identified 10,335 DEGs between the GL1 and GL2 libraries; 5,018 DEGs were upregulated and 5,317 DEGs were downregulated (Figs. 5B, 5C). Next, all DEGs were subjected to functional-enrichment GO analysis: 5,055 DEGs were significantly enriched in GO terms (Table S5) and were classified into three main categories, namely biological process, cellular component, and molecular function. In the biological process category, cellular process (2,715) made up the majority, followed by metabolic process (2,659) and single-organism metabolic process (2,421). In the cellular component category, organelle part (657) and intracellular organelle part (651) were predominant. In the molecular function category, catalytic activity (2,812) and binding (2,641) were prominently represented (Table S6).

We also performed an analysis of KEGG pathway enrichment in DEGs. Of the 10,335 DEGs, 1,188 unigenes were assigned and divided into 178 KEGG pathways (Tables S5, S7). Ribosome (123), lysosome (73), spliceosome (67), endocytosis (53) and plant-pathogen interaction (46) were highly represented pathways under salt stress (Table S7). Additionally, we identified the DEGs involved in the pathways of secondary metabolites biosynthesis and generated a heat map of these DEGs (Fig. 6). Phenylpropanoid biosynthesis (20) represented the largest proportion, followed by tropane piperidine and pyridine alkaloid biosynthesis (11), flavonoid biosynthesis (10), and isoquinoline alkaloid biosynthesis (eight). The annotation details of these DEGs are included in Table S8.

Genes involved in salt response

To identify genes responsible for salt response in Glehnia littoralis, we performed a BLAST search against the public databases. We found 1,661 unigenes encoding transcription factors (TF); 151 of these were DEGs (Table S5). These DEGs encoding transcription factors were divided into 36 types including MYB, basic leucine zipper (bZIP), ethylene responsive factor (ERF), helix-loop-helix (bHLH), and NAC. Among these, 71 DEGs were upregulated in response to NaCl treatment (Fig. 7). Additionally, unigenes of plant signal transduction pathways were collected from the DEGs according to the KEGG pathway. For example, 38 DEGs related to plant hormone signaling, 30 DEGs related to calcium signaling, and four DEGs related to phospholipase signaling were found to be responsive to NaCl treatment of Glehnia littoralis (Table S9). The heat map of these unigenes is displayed in Fig. 8.

Confirmation of DEGs by qRT-PCR

To confirm the accuracy and reliability of the transcriptome analysis data, 15 candidate unigenes associated with salt response and secondary metabolites were randomly selected for qRT-PCR assays (Fig. 9). These unigenes included predicted transcription factors and functional proteins. For example, MYB, bZIP, Whirkly, bHLH, and ERF transcription factors, calcineurin B-like protein (CBL), CBL-interacting protein kinase (CIPK), and phospholipases D and C were tested. Additionally, genes involved in second metabolites such as hydroxycinnamoyl transferase (HCT) and 4-coumarate-CoA ligase-like (4CL) were tested. As shown in Fig. 9, trends in the expression of these unigenes determined by qRT-PCR were consistent with those shown through RNA-Seq data, although the fold-changes detected by RNA-Seq and qRT-PCR did not match perfectly. Overall, these results suggest that our Glehnia littoralis transcriptome analysis using RNA-Seq was reliable.

Discussion

The medicinal plant Glehnia littoralis exhibits salt tolerance consistent with its original habitat. At present, little genomic information is available for Glehnia littoralis. In this study, we performed de novo transcriptome analysis of Glehnia littoralis using the Illumina HiSeq 2000 platform to investigate the salt tolerance mechanism and functional genes of Glehnia littoralis. Although only two libraries were constructed for de novo assembly and analysis with and without NaCl treatment, several plants were pooled for each sequencing treatment to reduce background difference among the samples as much as possible. Certainly, more sequencing replicates could be conducive to statistical screening and some outlier data could be removed by calculating the correlation among samples. Whatever, the expression results of the sequencing should be considered for validation by qRT-PCR and further gene function studies. In this study, our qRT-PCR results also supported the sequencing data. In total, we identified 10,335 DEGs and these DEGs were functionally categorised into a variety of physiological and molecular processes, which revealed the conserved mechanisms of salt responsive genes in Glehnia littoralis.

Plants subjected to salt stress display complex molecular responses such as stress gene expression, transcriptional regulation, and signal transduction networks (Deinlein et al., 2014). Therefore, we also focused on DEGs enriched in KEGG pathways under salt stress, and provided detailed information for genes in Glehnia littoralis associated with plant signal transduction pathways such as hormone signaling, calcium signaling, and phospholipase signalling (Fig. 8). ABA is an essential phytohormone that regulates plant stress response to environmental stimuli (Peleg & Blumwald, 2011). Drought and salinity cause osmotic stress in plants, and ABA is synthesised rapidly to trigger ABA-inducible gene expression, stomatal closure, and transpiration reduction to defend against water deficiency (Wilkinson & Davies, 2010; Yamaguchishinozaki & Shinozaki, 2006). In this study, we identified 524 unigenes involved in plant hormone signal transduction pathways. Of them, 34 DEGs including three pyrabactin resistance 1-like (PYL) ABA receptor genes were identified between the control and NaCl treatment samples (Fig. 8, Tables S7, S9). A previous study had shown that overexpression of PYL could enhance drought resistance and drought-induced leaf senescence in both Arabidopsis and rice by limiting transpirational water loss and water condition, thus generating an osmotic potential gradient leading water to preferentially flow to developing tissues (Zhao et al., 2016).

In addition, 30 DEGs were assigned to the calcium signaling pathway in Glehnia littoralis (Fig. 8, Tables S7, S9). Ca²⁺ signals are vital transducers and regulators in plant responses to environmental stimuli (Kudla, Batistic & Hashimoto, 2010). Typically, the salt-tolerance SOS3-SOS2-SOS1 pathway demonstrates the significance of Ca²⁺ signaling for excluding excess Na⁺ and maintaining cellular ion homeostasis (Qiu et al., 2002). SOS3 is a calcium sensor, belonging to the category of calcineurin B-like proteins (CBLs). SOS3 interacts with SOS2, a CBL-interacting protein kinase (CIPK), and the SOS3-SOS2 complex activates SOS1 (Na⁺/H⁺ antiporter) to promote the exclusion of Na⁺ to extracellular regions (Qiu et al., 2002). In Arabidopsis, AtCBL4 and CIPK24 correspond to SOS3 and SOS2, respectively (Gong et al., 2004). Furthermore, it has been demonstrated that AtSCABP8 (SOS3-like calcium binding protein8/CBL10) interacts with SOS2 and regulates SOS1 in the shoot response to salt stress (Quan et al., 2007). The plant CBL and CIPK gene families encompass several members; different CBL-CIPK complexes regulate different target proteins to respond to environmental stimuli and developmental needs (Kanwar et al., 2014; Liu et al., 2013; Pandey et al., 2015; Weinl & Kudla, 2009). In the current study, we assessed unigenes related to calcium signaling in G. littoralis. Some of these genes, such as CBLs and CIPKs, which are associated with plant salt-tolerance, were significantly upregulated. These data provide valuable genes targets for further research on physiological function in response to salt stress in G. littoralis. We also found that phospholipase-related DEGs were upregulated in G. littoralis under NaCl treatment (Fig. 8, Table S9). Phospholipases hydrolyse phospholipids. Phospholipase Ds (PLDs) and their hydrolysis product phosphatidic acid (PA), act as essential signal transducers and regulators involved in hyperosmotic stress responses (Hong et al., 2016; Zhang et al., 2014b). Previous studies have reported that several PLDs were activated by salt stress (Wang et al., 2014). In Arabidopsis, a lack of PLDa1 and PLDa3 resulted in more sensitivity to salt stress and a lower PA accumulation, and the change of PA level was suggested to affect plant salt response ( Hong et al., 2008; Yu et al., 2010). PA acts as a second messenger by its direct interaction with functional proteins, such as mitogen-activated protein kinase 6 (MPK6) and microtubule-associated protein (MAP65-1), to modulate the catalytic activity of binding proteins in response to salt stress (Yu et al., 2010; Zhang et al., 2012b). Moreover, a phosphoinositide-specific phospholipase C (OsPLC1)-mediated Ca²⁺ signaling pathway is essential for controlling Na⁺ accumulation in the rice response to salt stress (Li et al., 2017b).

Plant transcription factors respond to various environmental stimuli such as salt, drought, cold, and hormones and regulate gene expression by binding to core regions of promoters (Yadav et al., 2011). Several types of transcription factors associated with salt response have been identified in plants, including the MYB, NAC, bZIP, AP2/ERF, WRKY, and bHLH families (Cui et al., 2013; Hartmann et al., 2015; He et al., 2005; Jiang, Yang & Deyholos, 2009; Tao et al., 2011; Yang, Dai & Zhang, 2012; Zhang et al., 2009). In this study, we identified 151 differentially expressed transcription factors in G. littoralis under NaCl treatment (Table S5), suggesting their potential function in the regulation of G. littoralis salt response. For example, MYB constituted a large proportion of the DEGs encoding transcription factors in this study (Fig. 7). Previous studies have reported that the overexpression of the OsMYB3R-2, LcMYB1, and TaMYB56-B genes confers enhanced salt tolerance in transgenic plants (Cheng et al., 2013; Dai et al., 2007; Zhang et al., 2012a), and AtMYB20 and AtMYB73 act as negative regulators to enhance salt resistance (Cui et al., 2013; Kim et al., 2013). The changes in expression of these transcription factors in Glehnia littoralis could indicate their functional mechanism in salt stress.

G. littoralis roots (Radix Glehniae) are rich in polysaccharides, phospholipids, coumarins, coumarin glycosides and polyacetylenes (Tomshich et al., 1997; Yuan et al., 2002). G. littoralis grows naturally in coastal regions, and is considered a high-quality medicinal herb in traditional Chinese medicine. Previous studies have shown that salt stress could increase various secondary metabolites in plants, such as flavonoids, jasmonic acid, polyphenol, and polyamine (Ali, 2003; Ksouri et al., 2007; Mutlu & Bozcuk, 2007; Pedranzani et al., 2003; Ramakrishna & Ravishankar, 2011). Flavonoids were increased in response to salt stress in Hordeum vulgare seedlings (Ali, 2003). Ksouri et al. (2007) reported that polyphenol content and antioxidant activities were related to salinity in halophyte Cakile maritima. Free and bound polyamine content have been shown to be induced under salt stress in Helianthus annuus roots (Mutlu & Bozcuk, 2007). In tomato, endogenous jasmonic acid was increased under salt stress (Pedranzani et al., 2003). Based on the habitat of G. littoralis, we hypothesize that its active compounds may be associated with environmental salt stress. In a previous study of Peucedanum praeruptorum (Umbelliferae), transcriptome sequencing and high-performance liquid chromatography coupled with electrospray-ionization quadrupole time-of-flight mass spectrometry (HPLC-Q-TOF-MS/MS)-based metabolomics datasets were constructed for investigating the genes involved in coumarin biosynthesis and transport, as well for compound identification (Zhao et al., 2015). In total, 40,952 unigenes and 19 coumarin compounds were obtained and a few unigenes were predicted to be related to the formation of the coumarin core compounds in Peucedanum praeruptorum. Therefore, the G. littoralis transcriptome sequencing dataset could be used in future research to provide essential genes related to metabolism, and investigate biosynthetic pathways of the active compounds of G. littoralis.

Conclusions

In this study, we performed a comprehensive transcriptome analysis of G. littoralis in response to salt stress by Illumina 2000 sequencing. A significant number of unigenes, including DEGs, were identified and annotated. The DEGs dataset also provided candidate genes, which may be involved in salt tolerance or secondary metabolism pathways, to be used in subsequent functional analyses. These data will be helpful for future G. littoralis genomic studies, and will also be beneficial to studies on other species of Umbelliferae, and on halophytes.

Supplemental Information