De novo transcriptome sequencing of two cultivated jute species under salinity stress

Soil salinity, a major environmental stress, reduces agricultural productivity by restricting plant development and growth. Jute (Corchorus spp.), a commercially important bast fiber crop, includes two commercially cultivated species, Corchorus capsularis and Corchorus olitorius. We conducted high-throughput transcriptome sequencing of 24 C. capsularis and C. olitorius samples under salt stress and found 127 common differentially expressed genes (DEGs); additionally, 4489 and 492 common DEGs were identified in the root and leaf tissues, respectively, of both Corchorus species. Further, 32, 196, and 11 common differentially expressed transcription factors (DTFs) were detected in the leaf, root, or both tissues, respectively. Several Gene Ontology (GO) terms were enriched in NY and YY. A Kyoto Encyclopedia of Genes and Genomes analysis revealed numerous DEGs in both species. Abscisic acid and cytokinin signal pathways enriched respectively about 20 DEGs in leaves and roots of both NY and YY. The Ca2+, mitogen-activated protein kinase signaling and oxidative phosphorylation pathways were also found to be related to the plant response to salt stress, as evidenced by the DEGs in the roots of both species. These results provide insight into salt stress response mechanisms in plants as well as a basis for future breeding of salt-tolerant cultivars.


Introduction
Soil salinity is a major environmental stress imposed on plants that reduces agricultural productivity by restricting plant development and growth [1]. Salinity has primary effects including ion toxicity and osmotic stress as well as secondary effects such as oxidative stress [1]. Plants have a variety of salt tolerance mechanisms that depend on mitogen-activated protein kinase (MAPK/MPK) and hormone signaling as well as posttranslational modification of proteins.
Na + influx into roots occurs via different transporters. Plants use the Na+/H+ salt overly sensitive (SOS)1 antiporter, high-efficiency potassium transporter (HKT), and the tonoplastlocalized Na + , K + /H + exchanger (NHX) for sodium transport and detoxification [2,3] at 25˚C-28˚C in Yoshida nutrient solution in a greenhouse. And the salt tolerance of YY was only slightly better than that of NY. At the nine-leaf stage, six uniform seedlings of each species were selected; three were treated with 250 mM NaCl and the other three served as controls. After 12 h, the roots and leaves were collected for RNA extraction.

RNA isolation and transcriptome sequencing
Total RNA was extracted from tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, US) according to the manufacturer's protocol. RNA degradation and contamination was monitored in 1% agarose gels. RNA purity, integrity and the concentration was measured with a NanoPhotometer spectrophotometer (Implen, Westlake Village, CA, USA), RNA Nano 6000 Assay kit with the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA), and Qubit RNA Assay kit with the Qubit 2.0 flurometer (Life Technologies, Carlsbad, CA, USA). A total of 24 RNA sequencing libraries (NYCR, NY control root; NYCL, NY control leaf; NYSR, NY salt-stressed root; NYSL, NY salt-stressed leaf; YYCR, YY control root; YYCL, YY control leaf; YYSR, YY salt-stressed root; YYSL, YY salt-stressed leaf. three biological repetition for above each sample) were generated using the NEBNext Ultra RNA Library Prep kit for Illumina (New England Biolabs, Ipswich, MA, USA) according to the manufacturers' instructions. The quality of the libraries was assessed with the Agilent Bioanalyzer 2100 system. Libraries were sequenced on an Illumina HiSeq 4000 platform (San Diego, CA, USA) and paired-end reads were generated for transcriptome sequencing.

Transcriptome data analysis and annotation
Quality control of raw data was carried out with Perl scripts developed in house by removing reads containing adapter sequences, poly-N, and low-quality reads. Clean reads were assembled into a transcriptome using Trinity [17] with default settings for all parameters and mapped to transcripts; those with less than 5× coverage were removed. The transcripts which would be used as reference transcripts for analysis of DEGs were further assembled into non-redundant unigenes. Gene function was annotated using the following databases: National Center for Biotechnology Information (NCBI) non-redundant protein sequences (Nr); NCBI non-redundant nucleotide sequences (Nt); Protein Family (Pfam); Clusters of Orthologous Groups of Proteins (KOG); Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog (KO) database; and Gene Ontology (GO) using the Basic Local Alignment Search Tool with an E-value threshold of 10−5 [18] Analysis of differentially expressed genes (DEGs) Gene expression levels for each sample were estimated using RSEM [19]. Clean data from each sample were respectively mapped back onto the assembled transcripts, and the mapping results for each gene were used for differential expression analysis of unigenes. Differential expression analysis of the two conditions was carried out using the DESeq R package (1.10.1) [20]. Genes with P adjusted < 0.05 assigned by DESeq were considered as differentially expressed. GO enrichment and KEGG pathways analysis of DEGs were implemented with the GOseq R package [21] and KOBAS [22] software.

Quantitative real-time (qRT)-PCR analysis
Eight DEGs randomly selected from the RNA-seq results and eight DEGs selected from the list of some differentially expressed genes (DEGs) playing important roles in salinity stress were analyzed by qRT-PCR to validate RNA-seq results in the same samples used for RNA-seq with two independent biological and three technological replicates. The jute elongation factor-alpha (ELF) gene was used as the endogenous control [23]. Primer sequences for amplifying DEGs and ELF were listed in S1 File. For qRT-PCR, total RNA was reverse transcribed into firststrand cDNA using the M-MuLV Reverse Transcriptase kit (Fermentas, Burlington, ON, Canada) according to the manufacturer's instructions. The reaction was performed using the Fas-tStart Universal SYBR Green Master kit (Roche Diagnostics, Basel, Switzerland) according to the manufacturer's instructions in an optical 384-well plate on a GeneAmp PCR System 9700 (Applied Biosystems, Foster City, CA, USA). The cycle threshold method was used to calculate relative expression levels [24].  (Table 2). Transcripts were assembled into 72,278 non-redundant unigenes with an average length of 1240 bp. All sequencing data were deposited into the NCBI SRA database under accession numbers SRP116874 (included sequence data of 12 NY samples) and SRP119015 (included sequence data of 12 YY samples and unigenes annotation files). In total, 9108 (12.6%) unigenes were found to have homologs in the seven databases (Fig 1), with 52,417 (72.52%) annotated in at least one database. The largest number of unigenes (44,764,61.93%) had hits to protein sequences in Nr, while only 20,283 (28.06%) were mapped in KO.

Differentially expressed (D)TFs
In the study, a total of 2303 TFs were discovered (S6 File). There were more DTFs in YY than in NY (206/862 and 67/365 DTFs in leaves/roots of YY and NY, respectively) (Fig 2). In the  leaves, there were 32 DTFs common to the two species (S7 File and Fig 2). These DTFs belonged to 17 TF families, with homeobox (HB) and heat shock factor (HSF) being the most highly represented. In the roots, there were 196 DTFs (S8 File and
In order to identify pathways involved in the salt stress response, we carried out a KEGG pathways analysis of DEGs and found that pathways associated with metabolism and signaling were highly represented. There were 13 pathways enriched in both tissues in both species (Table 4) such as those related to plant hormone signal transduction and the peroxisome; and six pathways (Table 4) were enriched in roots in both NY and YY, including MAPK and Ca 2+ signaling. On the other hand, methane metabolism and cutin, suberin, and wax biosynthesis were only enriched in leaves of both species, although methane metabolism accounted for 10 and nine DEGs in NYL and YYL, respectively, and cutin, suberin, and wax biosynthesis accounted for four DEGs in the leaves of NY and YY (here, we only list the number of DEGs within enriched pathways in root and leaf of NY or YY in Table 4, while a few of DEGs (or no DEGs) without obviously enrichment in root and leaf of NY or YY were not showed (such as TCA cycle, only in root enrichment). Of the 127 DEGs common to both tissues and species, 13 were enriched in plant hormone signal transduction.

Some DEGs playing important roles in salinity stress
Here, numerous DEGs that played important roles in salinity stress were discovered. In Table 5, we list DEGs encoding HKT1, CBL-interacting protein kinases (CIPK), late embryogenesis abundant protein gene (LAE) etc. In the DEGs, the most of DEGs were the unigenes encoding CIPK (27), and 17 were up-regulated, following by the DEGs encoding MKK (9) with 6 down-regulated, 2 up-regulated and one up and down-regulated in NY or YY. The two DEGs (c122990_g1 and c122990_g2) were both up-regulated in the two tissues of both species. The two LAE DEGs and two SOS1 DEGs were all up-regulated in roots of YY. And one HKT1 was upregulated only in NY root.

Differential expression of plant hormone signal transduction genes
Hormones play an important role in the response to different environmental stressors (e.g., drought and salt stress) in plants. In this study, we identified numerous DEGs associated with plant hormone (ABA, auxin, cytokinin, gibberellin, ethylene, jasmonic acid, and salicylic acid) signaling in leaves and roots of both NY and YY, with ABA and cytokinin signaling being the most highly represented. In the ABA pathway, two genes encoding PYL were downregulated in all tissues and species whereas three other genes were downregulated in a subset of tissues and in one of the species. Five, one, and three genes encoding PP2C, SnRK2, and ABF were upregulated in all tissues and both species, respectively; and five genes encoding the three proteins were upregulated while one was downregulated in subset of tissues and in one of the species (Fig 4A). In the cytokinin signaling pathway, 13 genes encoding cytokinin response (CRE) 1 (5), B-APR (4), and A-type Arabidopsis response regulator (A-ARR) (4) were downregulated in at least one of four groups under salt treatment as compared to the control condition. Most genes encoding Arabidopsis histidine-containing phosphotransmitter (AHP) were upregulated in NYR, YYL, and YYR ( Fig 4B).

Pathways represented only in the roots of Corchorus
The Ca 2+ , MAPK signaling and oxidative phosphorylation pathways play major roles in the plant response to salt stress. Numerous DEGs identified in our study were associated with these pathways. Components of the Ca 2+ signaling pathway included calmodulin (CaM)/ CaM-like (CML) proteins and Ca 2+ -dependent protein kinase (CAMK) (Fig 5A). A total of 20 genes encoding CML (4/16 up-/downregulated) were found in YY or NY roots, with most exhibiting similar expression in the two species. One upregulated gene encoding CAMK was common to both YY and NY roots. In addition, nine DEGs encoding serine/threonine-protein phosphatase 2B catalytic subunit, four encoding Ca 2+ -transporting ATPase, and four encoding protein kinase (PK)A were detected in YY or NY roots. In the MAPK pathway, DEGs encoding HSP72 (n = 19), Ras-related C3 botulinum toxin substrate 1 (n = 12), MAPK1/3 (n = 9), and PP3C (n = 8) were detected in YY and NY roots (Fig 5B). In the oxidative phosphorylation pathway, 19 DEGs encoding H + -transporting ATPase subunits (V and F types) and six Common mechanisms of salt tolerance in Corchorus spp. encoding cytochrome c oxidase subunit were identified (Fig 5C). The expression profiles of DEGs in these three pathways were similar in the roots of the two species.

Validation of DEGs
In order to validate the RNA-seq results, we carried out qRT-PCR analysis for eight randomly selected DEGs and eight DEGs selected from Table 5 in YY and NY under salt stress. All of the DEGs showed differential expression under salt stress as compared to normal conditions, and the expression profiles were in accordance with the results obtained by RNA-seq (correlation coefficient = 0.85). Log2 (fold change) values for each DEG are shown in Fig 6.

Discussion
In this study, we carried out de novo transcriptome sequencing to investigate salt tolerance mechanisms in two jute species under salt stress. In previous studies, some reports showed that the number of DEGs was greater in roots than in leaves underlying salt stress. For example, under salt treatment, 4,884/3,692 up-/downregulated unigenes were discovered in cotton leaf samples; and much more up-/downregulated unigenes (6,303/11,068) were discovered in cotton roots DEGs [25]. while Yong et al. reported that [26] the number of DEGs was slightly greater in the roots (8,665 DEGs) than those in the leaves (7,795 DEGs) in Brassica napus. The number of DEGs in different tissues (e.g. roots and leaves) under salt stress might be related to crops. In the current study, much more DEGs common to both species in the roots than in leaves were discovered in jute, more similarly with in cotton; probably because that they were both classified as the Malvaceae [27]. Results of these studies might imply that roots were the main tissue responding to salt stress, especially in Malvaceae crops. In total, 127 DEGs were common to both tissues and both species, and the positive regulation DEGs (107) indwelled in all tissues and species were more than passive regulation DEGs (12) suggesting that the DEGs might play important roles in response to environmental stress in all of tissues and species. Most of the identified DEGs common to both tissues and both species were directly related to the response to environmental stress. For example, 13 DEGs were enriched in plant hormone signal transduction, and most of these were upregulated; DEGs involved in oxidationreduction and protein phosphorylation were also highly represented. Some DEGs, played important roles in salinity stress, were discovered. For example, there was evidence supports that HKT participated in the pathways of Na+ entry into the roots subjected to high salinity; and then Na+ was perceived and changed concentration of cytosolic Ca2+ decoded by Ca2 +-sensing proteins, such as CIPK. Oscillations of Ca2+ induced by salt stress were sensed by SOS3, a myristoylated Ca2+-binding protein; finally, the SOS1 was phosphorylated and activated through SOS pathway which extruded sodium in cells [28]. In addition, NHXs could compartmentalize Na+ in vacuoles to prevent toxic effects in the cytosol. In the study, numerous DEGs encoding HKT (1, up-regulated), CIPK (27, 17 up-regulated,), SOS1 (2, up-regulated) and NHX (3, up and down-regulated) were found [29]. The DEGs might be play crucial role in Na+ homeostasis. previous studies have shown that TF families such as AP2/EREBP, HB, MYB, NAC, bZIP, HSF, Orphans, FHA, and WRKY are linked to salt stress [30,31]. In the present study, DTFs expressed in leaves and roots differed; 17 DTF families including HB and HSF were detected in the former, while in the latter, 43 families including MYB, WRKY, and CCAAT were represented. Interestingly, 11 DTFs were common to both tissues in both species that are likely critical for the response to salt stress. Pathways obviously involved in the salt stress response were identified in the study. The pathways should be important for plants in response to abiotic stress. However, some Pathways were only enriched in roots or leaves in NY and YY, such as TCA cycle and Cutin, suberine and wax biosynthesis. For example, TCA cycle was obviously enriched in root of both species in this study, consistent with previous some studies. Kreuzwieser et al. proposed that hypoxia stress led to the inhibition of the TCA cycle and differential expression of some genes involved in TCA cycle only in roots of poplar trees [32]. And cutin, suberine and wax biosynthesis pathway was only enriched in leaves in NY and YY, in previous studies, the pathway was enriched in leaves [33] or in roots [34], but Lauter et al. reported cutin, suberine and wax biosynthesis term was enriched only in leaves of soybean under iron deficiency stress. suggesting in which tissues cutin, suberine and wax biosynthesis pathway was enriched might be influenced by numerous factors, such as species, stress conditions etc. Cutin, suberine and wax played important roles at plant environment interfaces by serving as a barrier controlling the movements of water and solutes. Which was thus important for the abilities of plants to withstand various abiotic stresses, such as drought and salinity [35]. So, which cutin, suberine and wax biosynthesis pathway was not enriched in root in NY and YY might result in uptaking more ion inside the root tissues from the rhizosphere and finally causes jute susceptibility under salinity stress.
The ABA signaling pathway is central to the salt stress responses in plants [1]. In this study, 20 DEGs were associated with this pathway, including genes encoding PP2C, SnRK2, and ABF (upregulated) and PYL (downregulated). Although PYL and PP2C expression is expected to be positively correlated based on what is known of ABA signaling, our observations are supported by recent studies [36,37]. The cytokinin signaling pathway includes cytokinin receptors; membrane-bound histidine kinases (e.g., CRE1) with a cytokinin-binding domain; AHP (which transmits signals to the nucleus); and B-and A-type ARR, which are phosphorylated by AHP in the nucleus [38,39]. Plants lacking or harboring mutant CRE1 show increased tolerance to high salt and drought conditions [40]; furthermore, B-and A-ARR mutations also increase salt tolerance [41], implying that cytokinin signaling negatively regulates the response to salt stress. In the present study, CRE1 and B-and A-ARR were downregulated in roots and leaves of both Corchorus species in response to high salinity. However, most AHP genes were upregulated, suggesting that they negatively regulate cytokinin signaling [42] to increase salt tolerance in plants.
CaM and CML are Ca 2+ sensors that transmit Ca 2+ signals to downstream effectors under conditions of environmental stress. Ca 2+ -loaded CaM/CML interacts with and regulates a broad spectrum of proteins including TFs, protein kinases and phosphatases, and metabolic enzymes [43]. Overexpression of CaM in Arabidopsis as a result of high salinity conferred salt stress tolerance via upregulation of the DNA-binding activity of MYB [44]. Interestingly, 20 DEGs corresponding to CML proteins were detected while MYB was the most highly represented TF family in YY and NY roots. CAMK upregulated in both YY and NY roots may enhance tolerance to salt stress [45] and interact with serine/threonine kinases that activate a Common mechanisms of salt tolerance in Corchorus spp. Na+/H+ antiporter, resulting in the extrusion of Na+ into soil or xylem for transporting to leaves [46]. CAMK and the serine/threonine protein kinase PP3C may activate Ca 2+ and MAPK signaling pathways in plants in response to high salt conditions or drought [10]. H + -ATPases were highly expressed in roots under salt stress, likely providing the proton-motive force required for maintaining membrane potential [47]. In the roots of both Corchorus species, there were 19 DEGs encoding H+-ATPase subunits under salt stress. This is consistent with previous studies [47] and demonstrates the importance of H + -ATPases for salt tolerance in plants.
In conclusion, our results provide a comprehensive resource for investigations into the mechanisms of salt response in plants, and can serve as a basis for breeding salt-tolerant cultivars of Corchorus species in the future.