Skip to main content
Download PDF

Transcriptome assembly, profiling and differential gene expression analysis of the halophyte Suaeda fruticosa provides insights into salt tolerance

BMC Genomics volume 16, Article number: 353 (2015) Cite this article

Abstract

Background

Improvement of crop production is needed to feed the growing world population as the amount and quality of agricultural land decreases and soil salinity increases. This has stimulated research on salt tolerance in plants. Most crops tolerate a limited amount of salt to survive and produce biomass, while halophytes (salt-tolerant plants) have the ability to grow with saline water utilizing specific biochemical mechanisms. However, little is known about the genes involved in salt tolerance. We have characterized the transcriptome of Suaeda fruticosa, a halophyte that has the ability to sequester salts in its leaves. Suaeda fruticosa is an annual shrub in the family Chenopodiaceae found in coastal and inland regions of Pakistan and Mediterranean shores. This plant is an obligate halophyte that grows optimally from 200–400 mM NaCl and can grow at up to 1000 mM NaCl. High throughput sequencing technology was performed to provide understanding of genes involved in the salt tolerance mechanism. De novo assembly of the transcriptome and analysis has allowed identification of differentially expressed and unique genes present in this non-conventional crop.

Results

Twelve sequencing libraries prepared from control (0 mM NaCl treated) and optimum (300 mM NaCl treated) plants were sequenced using Illumina Hiseq 2000 to investigate differential gene expression between shoots and roots of Suaeda fruticosa. The transcriptome was assembled de novo using Velvet and Oases k-45 and clustered using CDHIT-EST. There are 54,526 unigenes; among these 475 genes are downregulated and 44 are upregulated when samples from plants grown under optimal salt are compared with those grown without salt. BLAST analysis identified the differentially expressed genes, which were categorized in gene ontology terms and their pathways.

Conclusions

This work has identified potential genes involved in salt tolerance in Suaeda fruticosa, and has provided an outline of tools to use for de novo transcriptome analysis. The assemblies that were used provide coverage of a considerable proportion of the transcriptome, which allows analysis of differential gene expression and identification of genes that may be involved in salt tolerance. The transcriptome may serve as a reference sequence for study of other succulent halophytes.

Background

Salinity affects about 400 million hectares of land worldwide due to excessive irrigation and continues to increase in parallel with the population. Salinity in soil and water has caused substantial economic losses, including an estimated $230 million for the Indus Basin in Pakistan and $2 billion for the Colorado River basin in the U.S. [1,2]. An estimated total of 200 million hectares of new cropland is needed to feed the rapidly expanding population but only 93 million hectares are available for expansion and farming of traditional crops [1]. Attempts have been made with conventional crops to breed salt tolerance; however, these crops can only tolerate limited amounts of salt in their systems. The potential of halophytes, the natural flora of saline habitats, has been under-examined until recently and their utilization may allow production of useful crops on salty soils.

Suaeda fruticosa, a succulent shrub in the family Chenopodiaceae, is an obligate halophyte that grows optimally at 300 mM NaCl and has the adaptation to reduce sodium buildup for long term survival [2]. This perennial halophyte has a strong ability to accumulate and sequester Na+ and Cl− without the aid of salt glands, bladder or trichomes [3]. It is a good source of high quality edible oil [4], has potential for antiophthalmic, hypolipidaemic and hypoglycemic medicinal purposes [5], and has economic usage as forage for animals [6]. S. fruticosa also could help in bioremediation and reclamation of soils contaminated with toxic metals [7] and salinity [8]. Field studies showed that this plant can remove about 2646 kg of NaCl per hectare from the soil each year [9]. At optimum (300 mM NaCl) salt treatment of this species antioxidant enzymes trigger stress response through the activation of H2O2 − mediated Ca2+ uptake to maintain Na+ homeostasis at the cellular or tissue level [2]. Calcium ions, responsible for the overall signaling network of growth and development of the plant, are accumulated in the cell cytosol with the increase of Na+ [10]. At higher salinities, a significant reduction in growth is prevalent which might be due to the maximum threshold of the plant’s ability to adjust to specific ion toxicity and osmotic capability. Physiological data analysis has led to reports of ion accumulation, osmotic adjustments, maintenance of pressure potential and growth and production of glycinebetaine as part of a salt tolerance mechanism [11]. Previous studies of the impact of salinity on S. fruticosa have linked salt tolerance to its ability to uptake K+ in order to maintain a higher K+/Na+ ratio in the shoots. Higher sequestration of sodium and chloride in the shoot vacuoles together with the ability to synthesize osmoprotectants such as glycinebetaine has been suggested to maintain a favorable water potential gradient and protect cellular structures. Similar to Suaeda fruticosa, the majority of halophytes do not have glands or external bladders to modulate their tissue ion concentration therefore it has been seen to be a good model genus for the study of salt tolerance [12].

Next generation sequencing allows differential gene expression analysis of gene alleles and spliced transcripts, non-coding RNA and others, which will lead to identification of differentially expressed and/or unique genes. In this transcriptome paper, we report the identification of genes that are induced or repressed in plants grown under optimal salt conditions in comparison to low salt conditions. We generated a data set of transcript sequences from the roots and the shoots of Suaeda fruticosa. The genes were compared for differential expression under the indicated treatments using the assembled transcriptome, and common and tissue-specific patterns of transcriptomic responses were also analyzed. This first transcriptome study of Suaeda fruticosa expands our knowledge on global gene expression data for salt-accumulating halophytes that do not have external bladders.

Results and discussion

De novo transcriptome assembly and assessments of expressed sequenced tags

Experimental design

To prepare for the transcriptome assembly and analysis, total RNA was extracted from shoots and roots of Suaeda fruticosa. These include biological triplicates of cDNA libraries for S. fruticosa roots from plants grown without salt (R000), roots with 300 mM optimal salt (R300), shoots with no salt (S000) and shoots with 300 mM optimal salt (S300). Total mRNA was purified using oligo dT and transcribed into cDNA libraries using TruSeq RNA Sample Prep Kit for Illumina 100 bp paired-end sequencing.

Sequencing method and quality assessment of the reads

A total of 335.3 million reads of 100 bp were generated by Illumina Hi-seq platform (Sequence Read Archive Accession Number: SRX973396). The reads were filtered using Trimmomatic to remove adapters, FASTX toolkit and Sickle program to remove low quality reads and discard reads based upon the threshold of length. A total of 84.58% of the reads were trimmed and filtered totaling to 283,587,292 reads.

To normalize and assemble RNA-seq reads for de novo assembly, digital normalization was used for 283.6 million reads. K-mer hash of 21 with coverage of 30X was built from a set of reads to correct redundancy issues, variations in sequences, and potential errors among the reads. Since some reads with sequencing errors may escape quality score-based filtering steps, the reads with potential errors are flagged for removal from the dataset to improve the de novo assembly. Sequencing errors can affect the assembly algorithms so it is best to eliminate the reads that have non-uniform k-mer coverage. Reads that have non-uniform k-mer coverage create a problem with the assembly therefore it is necessary to normalize the reads to a certain threshold. This threshold represents the approximate minimum for de novo assembly to work optimally and efficiently. Digital normalization [13] was applied to the total of 283,587,292 paired end reads with k-mer size of 21 and k-mer coverage cutoff of 30X. The retained reads were normalized to 99,577,045 (Table 1) to remove overabundant reads, reduce noise of the sequenced sample and decrease the overall percentage of errors (Figure 1). The effect of digital normalization is to retain nearly all real k-mers while discarding the majority of erroneous and redundant k-mers. This step allows reducing the reads and obtaining a transcriptome assembly much faster than and superior to the assembly based on the full data set without affecting the quality of the assembly. Because the genes in the transcriptome have different levels of expression, k-mer distribution will not show any peak at any k value.

Table 1 Statistics of reads
Full size table
Figure 1
figure 1

Plot of total read pairs versus kept read pairs after digital normalization algorithms. The true k-mer counts are kept using digital normalization to reduce computational memory and correct redundancy.

Full size image

To assemble the Suaeda fruticosa transcriptome we utilized a genome-independent reconstruction approach. The strategy involved building a de Bruijn graph made of overlapping subsequences or k-mers using Velvet [14]. The overlapping bases allow building a graph of all the sequences that then traverse a path guided by read and paired-end coverage [15]. The path through the graph is reported as transcripts. To assemble the contigs into scaffolds, we used a de Bruijn graph software, Oases [16]. K-mer sizes from 35 to 99 were chosen to generate the assemblies. We assessed the quality of the assemblies based on the total number of transcripts, open reading frames predicted using Transdecoder and highest mapping percentage of the reads using GSNAP. The number of sequences, N50 values, mean length of the sequences and total base pair length for the contigs, scaffolds and unigenes were also determined (Table 1). Among the individual k-mer values, transcript numbers associated with different k-mer values vary from 1 to 450,588. K-mer length is related inversely to the number of transcripts generated. The highest N50 (computed by sorting the contigs from largest to smallest and then determining the minimum set whose sizes total 50% of the assembly) generated for the transcripts is 1,755 generated by a k-mer length of 59. Predicted open reading frames range from 1 to 108,112 bp with the highest number of ORFs being generated with a k-mer of 41. The highest N50 for open reading frames belongs to an assembly with a k-mer of 65 with 1,272 bp. Mapping coverage for the assemblies range from 30.39% to 72.91%. The highest percentage of reads mapping back to the assembly belongs to assemblies with k-mers 41 and 45 with 72.91% and 72.61% mapped. Assembly k-45 contains a higher percentage of proper pairs aligned with 59.56% and higher N50 compared to Assembly k-41, therefore it was chosen to be the assembly for the succeeding steps.

Assembly k-45 contained 296,776 contigs from Velvet with a N50 length of 1548 bp and mean size of 928 bp. We selected contigs that were greater than 200 bp in length. The contigs were assembled into scaffolds using Oases and yielded 273,824 contigs with an N50 length of 1669 bp and mean size of 1012 bp. The shortest scaffold is 152 bp and the longest one is 14,046 bp. Using CDHIT-EST, scaffold sequences were assembled into clusters and Transdecoder was used to predict open reading frames. The sequences were further screened for adapter and vector contamination using NCBI VecScreen and the UniVec database. Those sequences that are longer than 200 bp were kept. Using this pipeline, we obtained 54,526 high-quality unigenes with an N50 of 957 bp. The size range of the unigenes is between 200 to 6639 bp. There are 12,914 unigenes comprising 23.7% of the total that have lengths of more than 1000 bp. The mean size of the unigenes is 763 bp (Table 2).

Table 2 Statistics of sequence assembly
Full size table

Functional annotation, gene ontology assignments and analysis

The unigenes assembled were used as query for annotation using BlastX searches based on sequence homologies to the databases of the National Center for Biotechnology Information (NCBI) non-redundant (nr) protein database, RefSeq, SwissProt UniProt and the Kyoto Encyclopedia of Genes and Genomes (KEGG) using BLAST2GO. The summary of top hit distribution similar to Suaeda fruticosa unigenes is illustrated in Figure 2A. The species distribution with the lowest e-value matching the best sequence alignment result showed that the S. fruticosa transcriptome sequences have 8697 unigenes (16%) matching to Vitis vinifera (grapes), 3818 unigenes (7%) matching to Theobroma cacao (cacao tree) and 3127 unigenes (5.7%) similar to Beta vulgaris (beet). The closest halophyte species matching with Suaeda fruticosa is Populus trichocarpa (poplar tree) with 2327 unigenes (4.3% matching). For poplar only the initial analysis of the draft genome has been completed; additional mapping and sequencing is ongoing. Some of the halophytes mentioned in this paper do not have full annotation of genes submitted to NCBI database and some only contain partial transcriptome information. Figure 2B summarizes the data distribution summary from the sequences from the assembled transcriptome.

Figure 2
figure 2

Top hit distribution of matched unigenes among different species generated from BLASTX. (A). Data distribution summary from BLAST2GO shows blast hits, mapping results and annotated sequences (B).

Full size image

Gene names and GO terms were assigned to the transcripts based on homologies with an E-value threshold of 10-10. The data distribution summary for these sequences is shown in Figure 2B. Annotated sequences utilize assigned functional terms to query sequences from GO terms based on the gene ontology vocabulary. Mapped sequences are those with retrieved GO terms associated with the hits obtained after a BLAST search. The search produced 36,668 annotated sequences among 54,526 total transcripts, comprising 67.25% of total sequences. There are 8972 sequences comprising 16.45% of the total transcripts that did not surpass the annotation threshold and 6881 sequences or 12.6% had hits in the databases but lack functional information. A large proportion has no significant sequence alignment or hits in any of the databases, comprising 13,349 sequences or 24.5% of total transcripts which suggests that they may contain novel sequences or a high number of Suaeda fruticosa specific transcripts or transcript portions such as orphan untranslated regions.

Gene ontology encompasses a dynamic library for gene and protein roles in cells. This includes three main categories: Biological process, referring to the biological objective of the genes or gene products; Molecular function, defined by the biochemical activity of the genes or gene products; and Cellular components, referring to the place in the cell where the gene product is active [17]. Figure 3 illustrates the gene ontology annotation of the total assembled unigenes from the de novo assembled transcriptome of Suaeda fruticosa using BLAST2GO.

Figure 3
figure 3

Gene Ontology Summary of total assembled ESTs using BLAST2GO. Distribution of Gene Ontology Annotation of Suaeda fruticosa transcriptome. The results are summarized as follows: (A) Biological Process, (B). Cellular component (C) Molecular Function.

Full size image

In the biological process category, genes related to stress make up 1229 of the total unigenes annotated (Figure 3A). The other dominant subcategories were protein modification (933 unigenes), structural development (930 unigenes) and DNA metabolic process (923 unigenes). The following subcategories include unigenes involved in carbohydrate metabolism (798 unigenes), nucleobase-containing compound catabolic process (443 unigenes), organelle organization (348 unigenes), reproduction (513 unigenes), ribosome biogenesis (319 unigenes), signal transduction (594 unigenes), single organism development (457 unigenes), translation (346 unigenes), transmembrane transport (556 unigenes), lipid metabolism (513 unigenes), cofactor metabolism (314 unigenes), and cellular amino acid metabolism (850 unigenes). Figure 3B illustrates the cellular component category, which has a dominant subcategory of plastid (1374 unigenes), plasma membrane (1003 unigenes) and protein complex (941 unigenes). The molecular function category was comprised of protein coding genes involved in ion binding (3061 unigenes), oxidoreductase (961 unigenes), and those responsible for redox reactions of the cell and kinases (809 unigenes) (Figure 3C). These gene ontology annotations represent a profile for gene expression of Suaeda fruticosa suggesting that this species has diverse protein coding genes comprising its structural, regulatory, metabolic and stress response mechanisms.

Differential expression analysis

To acquire counts data for differential expression analysis, samples of different treatments (0 mM and 300 mM NaCl treatments) were mapped to the newly generated reference transcriptome using GSNAP (Genomic Short-Read Nucleotide Alignment Program) which utilizes computational methods to detect variants and splicing isoforms in short reads through merging and filtering position lists from a genomic index. It also detects short and long-distance splicing including interchromosomal splicing using probability models or a database of known splice sites [18]. Conversion of bam files into count data was performed using BamBam [19] to summarize the number of reads mapped to each annotated feature. Differential expression calls were made using the EdgeR package. Normalization is applied to the treatments and tissue types to provide accurate differential expression rather than individual quantification. The EdgeR package adjusted the analysis taking into account sequencing depths represented by library sizes. Variations between biological replicates were clustered closely using a multidimensional scaling (MDS) plot similar to that shown in Figure 4 to check for variations among replicates and samples.

Figure 4
figure 4

Multidimensional Scaling Plot for the sequencing libraries. Multidimensional Scaling Plot (MDS) is designed to indicate sample relationship similarity. Shoots and roots of 0 mM and 300 mM NaCl with their biological replicates are analyzed. (Key: S000A- shoots 0 mM replicate A, S000B- shoots 0 mM replicate B, S000C- shoots 0 mM replicate C, S300A- shoots 300 mM replicate A, S300B- shoots 300 mM replicate B, S300C- shoots 300 mM replicate C, R000A- roots 0 mM replicate A, R000B- roots 0 mM replicate B, R000C- roots 0 mM replicate C, R300A- roots 300 mM replicate A, R300B- roots 300 mM replicate B, R300C- roots 300 mM replicate C, .bam (bam files)).

Full size image

The replicates of treatments and their tissue types from transcriptome analysis were used to produce a multidimensional scaling plot (Figure 4), which allows us to see a spatial configuration of how similar or dissimilar the different treatments and biological replicates of S. fruticosa shoot and root samples are. The relationship of shoot treatments is more closely clustered together in comparison to root treatments. The tight clustering of the shoot data points means there are fewer variations among biological replicates in comparison to the root treatments. Root samples, however, have greater variations among the treatments and their biological replicates. This indicates that root tissues show less consistency with expression of genes among treatments. Common dispersions were then estimated on the distributions of reads across genes. Each gene gets an assignment of a unique dispersion estimate, which is to be compared to a common dispersion. The biological coefficients of variation versus the abundance were plotted (Figure 5). This specifies relative abundance of each gene variation between RNA samples and also measurement error estimated by the sequencing technology. From this sample, it shows a common dispersion of 0.37 and BCV of 61.09%. This means that common variation shows overall variability across the genome for this dataset and the common variation square root indicates high coefficient of biological variation. The genewise dispersions show a decrease at low average log counts per million. It indicates that at low expression level of genes or transcripts, the variability of gene abundance is high.

Figure 5
figure 5

Biological coefficient of variation plot. Genewise dispersion plot for twelve libraries is indicated. Estimation of genewise BCV allows observation of changes for genes that are consistent between biological replicates and giving less priority to those with inconsistent results. Generalized linear model is used to determine the evidence of significant difference of counts for a transcript or exon across conditions. The BH method is used in this dataset to control false discovery rate.

Full size image

The analyses were concentrated on genes that are significantly different in expression levels in the optimum salt transcriptome as compared to the low condition transcriptome. Genes whose adjusted p-values were less than 0.05 using the BH method were considered differentially expressed [20,21]. The BH method known also as FDR (false discovery rate) by Benjamini, Hochberg, and Yekutieli enables the user to control the false discovery rate, the expected proportion of false discoveries amongst the rejected hypotheses. The false discovery rate is a less stringent condition than the family-wise error rate, so these methods are more powerful than the others. RNA-seq gene expression for Suaeda fruticosa is visualized as an MA plot (log ratio versus abundance plot) in Figure 6. The red dots highlight transcripts that are differentially expressed among biological replicates and treatments. There are 475 genes that are downregulated and 44 genes are upregulated with a p-value <0.05 and false discovery rate <0.05. The results are consistent with the physiological data of Suaeda fruticosa [11] where at 0 mM NaCl treatment, more genes are downregulated in comparison to optimal growth of 300 mM NaCl.

Figure 6
figure 6

Differential expression genes plot. Plot of LogFCs against average count size, highlighting the differentially expressed genes in red. From the samples and the replicates, there are 475 genes identified to be downregulated and 44 genes that are upregulated with p-value of <0.05 and FDR rate of <0.05.

Full size image

Gene annotation and identification of differentially expressed genes

The differentially expressed genes were annotated using Blast2GO software against NCBI non-redundant protein database with a cut-off E-value of 10-10. Enrichment analysis was performed for the biological functions of the identified DEGs. Among 519 differentially expressed unigenes, 44 of them are upregulated upon salt treatment and 475 are downregulated. These genes were identified from BLAST nr, SwissProt and UniProt databases and assigned with Gene Ontology terms in biological process, molecular function and cellular component categories.

The top hit species distribution of these differentially expressed genes included grapes (Vitis vinifera) with 48 unigenes, orange (Citrus sinensis) with 35 unigenes, and Theobroma cacao with 29 genes. The closest halophyte is Populus trichocarpa with 13 unigenes and Mesembryanthemum crystallinum with 9 unigenes. Draft genome projects for both P. trichocarpa and M. crystallinum are currently ongoing while other halophytes only have partial transcriptome information available in the NCBI database.

From 519 differentially expressed genes, 391 unigenes have significant BLAST hits (75%) and the remaining 25% do not have any significant sequence alignments, which suggests that they might be genes that are novel or have not been reported in any other plant databases. There were 371 annotated sequences (71.5%), and 282 have InterProScan matches from the European Bioinformatics Institute (EBI) that can be mapped to GO terms and annotations while 162 of them have been assigned to gene ontology IDs. The summary in Figure 7 shows how many genes are assigned to at least one GO term and grouped into three main GO categories: biological processes (A), cellular component (B), and molecular function (C). Direct GO terms from Blast2GO were performed by counting annotated sequences in each term and suggesting the top terms. Among these sequences, 177 total unigenes are identified in the molecular function category of GO annotation. The top hits included genes functioning in ion binding (147 unigenes), kinase activity (43 unigenes) and DNA binding (40 unigenes). In the cellular component category, the top hits are genes found to be active in the nucleus (122 unigenes), protein complex component (121 unigenes) and plasma membrane (90 unigenes). For the biological process category, 205 unigenes have been assigned with GO terms and GO IDs. There are 174 unigenes that are important in biosynthetic process, 139 unigenes responding to stress and 124 unigenes involved in cellular nitrogen compound metabolic process. The differentially expressed genes were assigned to KEGG to identify pathways that these genes might be involved in related to salt tolerance. Among the annotated differentially expressed unigenes, the top hit included 6 sequences that are involved in both nitrogen and histidine metabolism. Others function in lysine degradation, glycerolipid metabolism and linoleic acid metabolism. Other pathways are illustrated in Additional file 1.

Figure 7
figure 7

Summary of differentially expressed ESTs using BLAST2GO. Differentially expressed transcripts were classified into 3 main GO annotations: Biological Processes (A), Cellular Component (B) and Molecular Functions (C). There are 25 GO terms for biological processes, 37 GO for molecular function and 15 GO for cellular component. A majority assigns the GO from biological process as stress response genes, genes responsible for oxidation-reduction and structure development. A few transcripts reflect oxidoreductase and kinase activity for molecular function. A majority of the transcripts is distributed to the nucleus and plasma membrane.

Full size image

Relative gene expression validation using qRTPCR analysis of RNA from 0 mM and 300 mM samples

To validate the results from the transcriptome analysis, we selected seven differentially expressed genes with putative functions related to salt tolerance. Specific primers were designed and optimized using PCR for the selected DE genes and for alpha tubulin as the endogenous control (Additional file 2). We amplified a cDNA library from six samples of 0 mM treated plants and six 300 mM treated samples. Analysis of transcript levels by qRTPCR showed that expression for all seven gene targets selected correspond with the differential expression patterns determined from the transcriptome analysis. Four targets (zeaxanthin epoxidase, aquaporin TIP2, dehydration responsive protein and glutathione S-transferase) show upregulation of mRNA expression while the other three targets (nitrate reductase, putative protein phosphatase and calcineurin B-like (CBL) 4–1) show downregulation upon 300 mM NaCl treatment compared to the absence of salt treatment (Figure 8).

Figure 8
figure 8

qRTPCR validation of the transcriptome data. Each panel shows the qRTPCR results for seven test genes. The annotated putative genes are listed on the x-axis and the mean fold change represented by the 2-ΔΔCT method relative to 0 mM treated samples are shown on the y axis. Error bars depict the standard error of the mean for 3 biological replicates. Significant differences (p < 0.05) are denoted with an asterisk and highly significant differences with p-value of <0.005 are represented with double asterisks.

Full size image

Putative salt tolerance-related genes

BLAST analysis data identified a large number of differentially expressed genes and we have grouped them in the following categories: 1. Genes responsible for enzymes, transcription factors, hormones, photosynthetic genes, detoxifiers and osmolytes for general metabolism, 2. Genes functioning as transporters for water and ion uptake, 3. Genes involved in regulation such as kinases and phosphatases, and 4. Genes that function to protect the cells against abiotic stress such as late embryogenesis abundant protein, heat shock proteins, osmoprotectants such as dehydrins and osmotins. The number of transcripts reported to be differentially regulated or expressed depends on the conditions being compared. In this study, we are comparing transcript expression between 0 and 300 mM NaCl treatment and their biological replicates. Upregulated genes are those with significant increased expression when treated with salt (300 mM). Those downregulated are annotated sequences with decreased expression with salt treatment. A summary of these sequences, their definitions and putative functions, and references from other halophytes or plants is shown in Table 3.

Table 3 Summary of selected differentially expressed genes in Suaeda fruticosa
Full size table

General metabolism genes

Genes that are involved in transcription, translation and post-translational modification have been seen to play roles in salt tolerance processes. WRKY transcription factors are important regulators for signaling mechanisms that modulate various plant processes. It has been found to interact with protein partners, MAP kinases, calmodulin, histone deacetylases, resistance proteins for autoregulation and transcriptional reprogramming [22]. It has also been suggested to be crucial for salinity tolerance [23]. From the differential expression analysis of the transcriptome, we have found WRKY transcription factor 72 to be significantly upregulated while WRKY transcription factor 6-like and WRKY DNA-binding protein isoform 2 are downregulated. The salt tolerant grass Festuca rubra ssp litoralis was found to have a differentially regulated WRKY-type transcription factor in response to salinity [24]. Transient expression studies have also found OsWRKY72 and OsWRKY77 to be activators of ABA signaling but repressors of gibberelic acid signaling in aleurone cells [25]. Moreover, in Arabidopsis, AtWRKY6 negatively autoregulates its own promoter to influence senescence and pathogen defense-associated PR1 promoter activity. This targets SIRK, a gene encoding a receptor-like protein kinase that is strongly induced during leaf senescence. The activation of SIRK is dependent on WRKY6 function [26]. These studies suggest that WRKY72 trancription factor is upregulated to respond to ABA signaling, important for stress tolerance while downregulating protein-coding genes involved in senescence for protection and defense. Gibberelic acid (GA) genes, which regulate many aspects of growth and development of plants, are involved in the synthesis of gibberellin hormone. In Arabidopsis, reduction of bioactive GA is shown via an increase in gibberellin 2-oxidase 7 (GA2ox7). This leads to accumulation of DELLA proteins, which are transcriptional regulators that repress GA-responsive growth and development, inhibiting plant growth [27]. Downregulation of GA2ox2 is observed at 300 mM salt treatment in Suaeda fruticosa. This suggests that the decrease deactivates bioactive GA [28]. GA genes were regulated at 200 mM NaCl in S. europaea similar to homologues of gibberellin 3-oxidase and gibberellin 20-oxidase in P. trichocarpa. Two DELLA domain GRAS family transcription factors were downregulated in plants treated with 200 mM salt [29].

Both 40S ribosomal protein S4 and 60S ribosomal protein L18-2-like that are upregulated in S. fruticosa are part of a group of SRP-dependent co-translational proteins targeting to membranes responsible for translation and protein binding [30]. Ribosomal protein 40S and 60S and RNA binding family protein are also highly upregulated in this transcriptome study. Similar studies were performed in Poplar euphratica, which found that ribosomal 60S rRNA, important for plastidic and nuclear protein synthesis, is increased in response to salinity [31]. 60S acidic ribosomal protein P2, known to play an important role in the elongation step of protein synthesis and other RNA-binding family proteins are upregulated in 300 mM NaCl treated Suaeda fruticosa. However, the gene encoding pre-mRNA processing protein 40C undergoes downregulation in salt treated plants. This protein has been found to be a coactivator involved in regulated transcription of RNA polymerase II-dependent genes important in transcription and other regulatory mechanisms [30]. Some DNA binding proteins also show concerted regulation upon salt treatment. DNA-binding escarola-like protein responsible for late flowering and leaf development and F-box kelch repeat protein AT1g80440-like are upregulated while MADS-box transcription factor AGL24, an early target of transcriptional repression at floral transitional stage, DNA-binding protein with zinc finger isoform1 and F-box protein AT1g78280-like transferase involved in regulation of transcription are downregulated.

An increase in reactive oxygen species causing damage to cellular components is evident when salinity increases. Genes that are responsible in regulating redox reactions are usually involved in protecting the cell environment during these stresses [32]. Upregulation of glutathione-S-transferase tau 1 (GST) and glutathione transferase were seen in the differential expression analysis of S. fruticosa. Similarly, glutathione S-transferases were greatly increased upon salt treatment in roots of the halophyte Salicornia europaea [33], Suaeda maritima and Reaumuria trigyna [32,34]. The Suaeda salsa GST gene was introduced into Arabidopsis and improved salt tolerance after overexpression in transgenic plants. Glutathione content increased in salt-stressed Arabidopsis and promoted a higher level of salt tolerance [35]. The level of glutathione is increased at 0 mM and 900 mM treatment and decreased at the optimal condition of 300 mM NaCl in S. fruticosa [2].

A similar trend of higher salt tolerance is seen in tobacco seedlings upon overexpression of GST and these genes have been found to be responsible for increased protection against toxins [36]. Some proteins important for seed production and growth show differential expression in S. fruticosa. Germin-like protein, found to be an important plant marker for salt stress regulation and suggested to undergo change when salt-tolerant plants are subjected to salt stress has been found to be significantly upregulated upon salt treatment [37,38]. An ortholog of flowering promoting factor 1-like protein 3, which promotes flowering in Arabidopsis, and auxin-induced protein 5NG4-like gene involved in transport of molecules functioning downstream of the auxin response and responsible for root formation are also upregulated [39]. Some genes encoding proteins involved in protection such as pathogenesis-related protein, chitinase, peroxisomal ascorbate peroxidase (APX), and plant cadmium resistance 2-like are also increased. Plant chitinase plays an important role in plant defense and enhances resistance and tolerance to heat, salt and drought [30]. Overexpression of chitinases in transgenic tobacco has been shown to enhance biotic and abiotic stress tolerance [40]. In tobacco cells, APX functions in the regeneration of NAD+ and is usually induced by high temperature stress and functions against toxic reactive oxygen species [41]. In the halophyte Atriplex halimus L., chloride salinity reduces cadmium accumulation as salinity resistance is found to be closely associated with the gene loci responsible for cadmium extraction [42-44]. Proteins containing chaperone domains and DNAJ-16 like chaperon protein are also decreased upon salt treatment. The DNAJ protein family is included in the group of heat shock proteins functioning as molecular chaperones, and is associated with HSP70 and involved in resisting environmental stresses in Suaeda maritima [34]. Specifically DNA-J16 in Arabidopsis is encoded by the gene known as Altered Response to Gravity 1 (ARG1), and mediates gravity signal transduction and hypocotyl gravitropism [45,46]. Other genes that are downregulated include ethylene-responsive transcription factor rap2–7 like isoform, senescence-associated protein, stem specific protein TSJT1-like, root hair protein F3H11–7, cell wall protein AWA1-like isoform X1 for cell wall organization, and callose synthase 7, a major component of pollen tubes and pollen cell walls. Molecular mechanisms of cellular calcium changes have been seen with the downregulation of calcineurin B-like protein (CBL) and calmodulin binding isoform 1 upon salt treatment suggesting their potential role as regulators of salt and drought responses [47]. Calmodulin mediates auxin signaling and responds to stresses in Arabidopsis [48]. CBL interacts with CIPK serine-threonine protein kinases and mediates activation of AKT1 in response to low potassium conditions and stomatal movement [49].

Various photosynthetic genes have been found to be differentially upregulated upon salt treatment in S. fruticosa. These include genes encoding photosystem II protein z, d2 protein, cp43 chlorophyll apoprotein, chloroplast RF19, zeaxanthin epoxidase, chloroplastic like isoform X2 and sufE-like chloroplastic protein. Significant induction has also been found in the halophyte Salicornia europaea in which photosynthetic genes, PSI and PSII pigment binding proteins, b6f complex and ATPase synthase CF1 are upregulated in salt treated plants [33]. Some genes encoding light-induced chloroplastic protein, CRS2-associated factor, thioredoxin-like protein chloroplastic like, triose phosphate chloroplastic-like isoform X2, probable chlorophyll b reductase chloroplastic-like and phosphate chloroplastic-like are downregulated in S. fruticosa. While some of these proteins have no definite functions determined yet, chlorophyll b reductase has been found to play a role in maturation and storability of seeds in Arabidopsis. Arabidopsis plants lacking chlorophyll b show a stay-green phenotype in leaves [50]. This suggests that as chlorophyll b reductase decreases in plants, they tend to prevent chlorophyll degradation.

Ion transporters (transporters and aquaporins)

Homeostasis of the cellular environment involves the maintenance of cellular uptake to control ionic balance. Since a large influx of extracellular Na + occurs in halophytes, plants require high amounts of K+ (100–200 mM) to lower the amount of Na + and maintain osmosis [51]. Aquaporin tonoplast intrinsic proteins showed upregulation in salt treated S. fruticosa. Aquaporins are membrane proteins that facilitate uptake of soil water and mediate regulation of root hydraulic conductivity. They are also involved in compartmentalization of water and are found in halophytes to play a role in maintaining osmosis and turgor of plant cells [52]. The halophyte Schrenkiella parvulla contains high numbers of aquaporins for tolerance to boron toxicity [53]. In Poplar species, some aquaporins are decreased to prevent water loss during salt stress [31]. Some other transporters that are upregulated upon salt treatment include high-affinity nitrate transporter 3.1-like and nitrate transporter 1.5 important for nitrate uptake, aluminum-activated malate transporter 10 for increased aluminum tolerance, flavonol 4-sulfotransferase for auxin transport, bidirectional sugar transporters and glucosyltransferase for glucose and other sugar transport, seed storage/lipid transfer protein responsible for metabolism and transport, and ATPase subunit 1. Other halophytes such as Schrenkiella parvula and Thellungiella showed upregulation of genes encoding for ATPases that are necessary for large influx of ions [53,54].

Studies have shown vacuolar Na+/H+ antiporter to be important for salt tolerance through Na + sequestration [55]. However, in Suaeda fruticosa, sodium transporter HKT1-like is shown to be downregulated. In Arabidopsis, HKT1 knockouts accumulate the highest concentration of Na + in the shoots suggesting a role in maintenance of Na + concentration [56]. Some other ion transporters are also downregulated such as sodium pyruvate chloroplastic co-transporter, magnesium transporter NIPA2, vacuolar iron transporter, vacuolar proton ATPase A1-like and ASNA1 (arsenic pump driving ATPase). Glutamyl-tRNA amidotransferase involved in carbon-nitrogen ligase activity, tonoplast dicarboxylate transporter-like protein for malate transmembrane transport and regulation of intracellular pH and galacturonosyltransferase 12-like for glycan and pectin biosynthesis are also decreased with salt treatment.

Regulatory molecules (kinases and phosphatases)

Differentially regulated molecules such as kinases and phosphatase are involved in regulation of proteins involved in osmolyte synthesis and detoxification by oxidants. They are suggested to play a role in ionic and osmotic homeostasis and modulate ion transport for salt tolerance [57]. Cysteine-rich receptor like protein kinase, phosphatase 2C family protein including phosphatase 2C 15-like isoform X1 and purple acid phosphatase 27-like are upregulated at 300 mM NaCl treatment. Protein phosphatase 2C (PP2C) regulates signal transduction pathways. In Thellungiella, A-type PP2C phosphatases are generally upregulated in response to abscisic acid [58]. Moreover, there are other kinases that are downregulated in this study such as CDPK-related kinase 1, PERK1 kinases, casein kinase I2-like protein, and serine-threonine protein kinase HT1 and phosphoinositide 4-kinase gamma 4. Serine threonine protein kinase HT1 is important for regulation of stomatal movement in response to carbon dioxide [59] while CDPK-kinase 1 has been shown to play an important role in mediating signal transduction of growth and development [30]. In rice, OsCDPK1 negatively regulates the expression of enzymes for gibberellic acid biosynthesis. This also transduces post-germination of Ca2+ signal from sugar starvation and gibberellic acid to prevent drought stress injury [60]. Some phosphatases are also downregulated such as serine-threonine protein phosphatase pp1-like, phosphatase 2C 76 isoform 1 and PP2A catalytic subunit. In Arabidopsis, transcription factor MYB20 negatively regulates 2C serine-threonine protein phosphatases to enhance salt tolerance [61].

LEA genes

Late embryogenesis abundant (LEA) proteins comprise a group of proteins that have crucial roles in cellular dehydration tolerance. They have been associated with tolerance to dehydration caused by freezing, salinity or drying. During stress conditions such as salinity, plant hormone abscisic acid (ABA) is produced to develop tolerance against drought. Some genes are induced to trigger the production of ABA [62].

Overexpression of LEA proteins can improve stress tolerance of transgenic plants [62]. In this transcriptome study, salt treatment causes upregulation of dehydration-responsive RD22-like protein. RD22 expression in Arabidopsis is mediated by abscisic acid (ABA). This is also induced by salt stress and dehydration [63] and is expressed during early and middle stages of seed development. Housekeeping gene HSP20 chaperone superfamily is found to be downregulated upon salt treatment. HSP20 family has been associated with the most stress-general expression pattern including salt stress in Arabidopsis [64].

Conclusions

This study provides an overview of the genes present in a non-model plant species and identifies the genes associated with salt tolerance. The assembled transcriptome was used for differential expression studies and gene annotations. We have identified 519 genes that are differentially expressed based on p-value and adjusted false discovery rate of less than 0.05. The same pattern of differential expression for seven of these genes was confirmed by qRT-PCR analysis, which each showed similar levels of up- or down-regulation (Figure 8).

The annotation of genes using next generation sequencing is more readily available through the advancement of technology. Analysis of predicted genes allows assumptions to be made on the complexity of genetic mechanisms for this plant. RNA sequencing generates an enormous amount of data in terms of identifying the transcripts, however the challenges remain with the analysis. One of the major problems is the development of expression metrics that will allow comparisons of different expression levels and also provide identification of differentially expressed genes. We have utilized a combination of approaches to conduct this analysis for the Suaeda fruticosa transcriptome. The reference transcriptome assembly was not previously available and this species does not have any close relative plant that can serve as a basis for the expression analysis.

This study reports comprehensive information about the transcriptome of the succulent halophyte S. fruticosa. This will provide a basis for further study of the mechanism of salt tolerance, discovery of novel genes involved and comparison of expression profiles with no salt and optimal salt concentration. The de novo transcriptome generated in this study provides a useful source of reference sequence for succulent halophytes.

Methods

Plant materials and RNA isolation

Seeds of Suaeda fruticosa obtained from the Institute of Sustainable Halophyte Utilization, University of Karachi, Pakistan were planted and grown at Brigham Young University, Provo, Utah, U.S.A. according to protocol [2]. Plant samples of 100 mg of frozen plant tissue from roots and shoots of low (0 mM NaCl) and optimal (300 mM NaCl) salt conditions were ground in liquid nitrogen to a fine powder. Total RNA was extracted from these tissues using a Trizol-based method or QIAGEN RNeasy Mini kit. The RNA was analyzed for quality and concentration using the Agilent Technologies 2100 Bioanalyzer. High quality total RNA samples should give two distinct peaks and yield an RNA Integrity Number (RIN) value greater than 8.

Illumina sequencing platform

The Illumina RNA-Seq library preparation protocol includes poly-A RNA isolation, RNA fragmentation, reverse transcription to cDNA using random primers, adapter ligation, size-selection from a gel and PCR enrichment [65]. The batch of libraries was sequenced at the BYU sequencing center and by Otogenetics (Norcross, GA) using Illumina Hi-seq 2000 sequencer. This includes cDNA libraries of Suaeda 0 mM NaCl-treated shoot and roots in triplicates and cDNA libraries of Suaeda 300 mM NaCl-treated shoot and roots in triplicates. The paired-end library was developed according to the protocol of the Paired-End sample Preparation kit (Illumina, USA).

Bioinformatics analysis

Quality trimming and digital normalization

The adapters of raw RNA sequence data were trimmed using Trimmomatic v.0.27. FastX toolkit and Sickle paired-end trimmers were used to determine low quality reads towards the 3’ and 5’ ends of the reads. The software for digital normalization is available electronically through http://ged.msu.edu/papers/2012-diginorm/. A python script was used to interleave the paired-end reads file https://github.com/ged-lab/khmer/tree/2012-paper-diginorm/sandbox. Khmer software package available at https://github.com/ged-lab/khmer/tree/2012-paper-diginorm was used to perform three-pass normalization steps. Loading sequences needed for khmer software works with screed packages through https://github.com/ged-lab/screed/ (khmer and screed are ©2010 Michigan State University, and are free software available for distribution, modification, or redistribution under the BSD license). The details of quality trimming and digital normalization are available in Additional file 3.

De novo assembly and gene ontology

The high quality concatenated reads of shoots and roots (Sequence Read Archive Accession Number: SRX973396) were assembled using Velvet v. 1.2.10 and Oases v. 0.2.08 with optimized determined k-mer size of 45 with an average insert length of 300 bp and minimum contig length of 200. Further screening of sequences for vector contaminations and adapters were performed using NCBI Univec and VecScreen. Only assembled transcripts longer than 200 bp were kept. De novo assembly scripts are available in Additional file 3. We ran the clustering methods using CDHIT-EST v.4.5.4–2011–03–07 on the assembly. All Illumina assembled unigenes were searched against nr database in NCBI, Swiss-Prot, UniProt, and Kyoto Encyclopedia of Genes and Genome (KEGG) with the BLASTX algorithm. The E-value cut-off was set to 10-10. Genes were identified according to best hits against known sequences. Prediction of gene ontology (GO) terms, sequences functions, metabolic pathways in KEGG databases were performed.

Sequence analysis

The assemblies were transferred into Transdecoder, an open-reading frame predictor software under the Trinity package, which reports candidate coding regions within the transcripts. For each assembly, the number of transcripts, N50 values and the total length of the assemblies are identified. The analysis of the efficiency of assemblies is performed using GMAP and GSNAP v. 2013–11–27. GMAP maps and aligns cDNA sequences originally used for genomic mapping then GSNAP aligns single-end or paired-end reads. It can detect short and long distance splicing using probabilistic models or database of known splice sites.

Differential expression analysis

To determine the DEGs (differentially expressed genes) between different treatments of shoots and roots of Suaeda fruticosa, gene expression level analysis was performed using the EdgeR package from R [66]. Calculated gene expression can be directly used for comparing the differences in gene counts between treatments and tissue types. Generalized Linear Models were used for data analysis to take account of different salt conditions and tissue types of biological replicates. This determines the evidence of significant difference of counts for a transcript or exon across experimental conditions. The estimation for biological variation is measured. DEGs were identified and subject to further annotation using BLAST2Go.

Validation of differentially expressed genes through qRTPCR

Several putative annotated genes were selected for validation of differential expression using qRTPCR. These include aquaporin TIP2, protein phosphatase, calcineurin b-like protein (CBL) 4–1, zeaxanthin epoxidase, dehydration responsive protein, glutathione S-transferase and nitrate reductase. We selected alpha tubulin as an endogenous control. The primers for these genes were designed from the Suaeda fruticosa transcriptome sequences and optimized for PCR (Additional file 2).

For each qRTPCR reaction, 1 ug of RNA of 0 mM and 300 mM NaCl treated samples were reverse transcribed into cDNA using oligo (dT) primers, and the cDNA libraries produced were used for qRTPCR using the method of Haddad et al. [67]. To assess validation for each gene, qRTPCR data were analyzed based on ΔΔCT and 2-ΔΔCT method [68]. The ΔCT value of each gene was calculated by subtracting the CT value of the endogenous control from the CT value of the target gene. Each gene’s mean ΔΔCT value.

2-ΔΔCT and standard error of the mean were calculated using the data analysis package in Microsoft Excel. Data were plotted as mean fold change (2-ΔΔCT). Significant differences (p < 0.05) were determined using a one-tailed two sample t-test assuming equal variances for comparison of the fold change values between groups using GraphPad software.

Availability of supporting data

Raw Illumina sequences are available at NCBI Sequence Read Archive under Suaeda fruticosa accession SRX973396 (SRA run accessions: SRR1946790, SRR1946833, SRR1946834, SRR1947647, SRR1947648, SRR1947649, SRR1947655, SRR1947656, SRR1947661, SRR1947683, SRR1947686, SRR1947688). Transcriptome sequence information is deposited in the Transcriptome Shotgun Assembly Sequence Database: BioProject ID: PRJNA279962 and PRJNA279890. The following supplementary information files are publicly available at LabArchives (DOI: http://0-dx-doi-org.brum.beds.ac.uk/10.6070/H4Z03650). Additional file 1 shows the exported name list of assembled sequences, their annotated names, gene ontology terms, associated enzymes and KEGG maps generated by BLAST2GO. Additional file 2 shows the summary of the number of unigenes assigned to different pathways. Additional file 3 lists the codes used for Bioinformatics analysis. Additional file 4 lists the primers used for qRTPCR on selected genes.

Abbreviations

R000:

Roots treated without salt

R300:

Roots treated with 300mM NaCl

S000:

Shoots treated without salt

S300:

Shoots treated with 300mM salt

ORF:

Open reading frame

nr:

Non-redundant database

GO:

Gene ontology

KEGG:

Kyoto Encyclopedia of genes and genomes

BLAST:

Basic local alignment search tool

GSNAP:

Genomic short read nucleotide alignment program

MDS:

Multidimensional scaling plot

BCV:

Biological coefficient of variation

FDR:

False discovery rate

BH:

Benjamini, Hochberg

GA:

Gibberelic acid

GST:

Glutathione S-transferase tau 1

APX:

Peroxisomal ascorbate peroxidase

PP2C:

Protein phosphate 2C

LEA:

Late embryogenesis abundant proteins

ABA:

Abscisic acid

DEG:

Differentially expressed genes

References

  1. Beltran JM, Manzur CL. Overview of salinity problems in the world and FAO strategies to address the problem. Proceedings of the International Salinity Forum. 2005:311–313.

  2. Hameed A, Hussain T, Gulzar S, Aziz I, Gul B, Khan MA: Salt tolerance of a cash crop halophyte Suaeda fruticosa: biochemical responses to salt and exogenous chemical treatments. Acta Physiologiae Plantarum. 2012;34:2331–2340.

  3. Labidi N, Ammari M, Mssedi D, Benzerti M, Snoussi S, Abdelly C. Salt excretion in Suaeda fruticosa. Acta Biol Hung. 2010;61:299–312.

    Article  PubMed  Google Scholar 

  4. Weber DJ. Adaptive mechanisms of halophytes in desert regions. Salinity and water stress. 2008;44:179–85.

    Article  Google Scholar 

  5. Bennani-Kabachi N, El-Bouayadi F, Kehel L, Fdhil H, Marquie G. Effect of Suaeda fruticosa aqueous extract in the hypercholesterolaemic and insulin-resistant sand rat. Therapie. 1999;54:725–30.

    Google Scholar 

  6. Towhidi A, Saberifar T, Dirandeh E. Nutritive value of some herbage for dromedary camels in the central arid zone of Iran. Tropical Animal Health Pro. 2011;43:617–22.

    Article  Google Scholar 

  7. Bareen F, Tahira SA. Metal accumulation potential of wild plants in tannery effluent contaminated soil of Kasur, Pakistan: field trials for toxic metal cleanup using Suaeda fruticosa. J Hazard Mater. 2011;186:443–50.

    Article  CAS  Google Scholar 

  8. Khan MA, Ansari R, Ali H, Gul B, Nielsen BL. Panicum turgidum: a sustainable feed alternative for cattle in saline areas. Agric Eco Env. 2009;129:542–6.

    Article  Google Scholar 

  9. Chaudhri II, Shah BH, Naqvi N, Mallick IA. Investigations on the role of Suaeda fruticosa Forsk in the reclamation of saline and alkaline soils in West Pakistan plains. Plant Soil. 1964;21:1–7.

    Article  Google Scholar 

  10. Anil VS, Rahjkumar P, Kumar P, Mathew MK. A plant Ca2+ pump, ACA2, relieves salt hypersensitivity in yeast. Modulation of cytosolic calcium signature and activation of adaptive Na + homeostasis. J Biol Chem. 2008;283:3497–506.

    Article  CAS  PubMed  Google Scholar 

  11. Khan MA, Ungar IA, Showalter AM. The effect of salinity on the growth, water status, and ion content of a leaf succulent perennial halophyte, Suaeda fruticosa (L.) Forssk. J Arid Environ. 2000;45:73–84.

    Article  Google Scholar 

  12. Flowers TJ, Colmer TD. Salinity tolerance in halophytes. New Phytol. 2008;179:945–63.

    Article  CAS  PubMed  Google Scholar 

  13. Brown CT, Howe A, Zhang Q, Pyrkosz A, Brom T. A reference-free algorithm for computational normalization of shotgun sequencing data. arXiv: 12034802. 2012;

  14. Zerbino DR. Using the Velvet de novo assembler for short-read sequencing technologies. Curr Protoc Bioinformatics. 2010; Chapter 11: Unit 11 15.

  15. Garber M, Grabherr MG, Guttman M, Trapnell C. Computational methods for transcriptome annotation and quantification using RNA-seq. Nat Methods. 2011;8:469–77.

    Article  CAS  PubMed  Google Scholar 

  16. Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics. 2012;28:1086–92.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  17. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  18. Wu TD, Nacu S. Fast and SNP-tolerant detection of complex variants and splicing in short reads. Bioinformatics. 2010;26:873–81.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  19. Page JT, Liechty Z, Huynh M, Udall JA: BamBam: genome sequence analysis tools for biologists. BMC Research Notes. 2014;7:829.

  20. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practival and powerful approach to multiple testing. J Royal Statistical Soc Series. 1995;57:289–300.

    Google Scholar 

  21. Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Statistics. 2001;29:1165–88.

    Article  Google Scholar 

  22. Rushton PJ, Somssich IE, Ringler P, Shen QJ. WRKY transcription factors. Trends Plant Sci. 2010;15:247–58.

    Article  CAS  PubMed  Google Scholar 

  23. Garg R, Verma M, Agrawal S, Shankar R, Majee M, Jain M. Deep transcriptome sequencing of wild halophyte rice, Porteresia coarctata, provides novel insights into the salinity and submergence tolerance factors. DNA Res. 2014;21:69–84.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Diédhiou CJ, Popova OV, Golldack D. Comparison of salt-responsive gene regulation in rice and in the salt-tolerant Festuca rubra ssp. litoralis. Plant Signal Behav. 2009;4:533–5.

    Article  PubMed Central  PubMed  Google Scholar 

  25. Xie Z, Zhang Z-L, Zou X, Huang J, Ruas P, Thompson D, et al. Annotations and functional analyses of the rice WRKY gene superfamily reveal positive and negative regulators of abscisic acid signaling in aleurone cells. Plant Physiol. 2005;137:176–89.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  26. Robatzek S, Somssich IE. Targets of AtWRKY6 regulation during plant senescence and pathogen defense. Genes Dev. 2002;16:1139–49.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  27. Achard P, Gong F, Cheminant S, Alioua M, Hedden P, Genschik P. The cold-inducible CBF1 factor–dependent signaling pathway modulates the accumulation of the growth-repressing DELLA proteins via its effect on Gibberellin metabolism. The Plant Cell Online. 2008;20:2117–29.

    Article  CAS  Google Scholar 

  28. Alvey L, Harberd NP. DELLA proteins: integrators of multiple plant growth regulatory inputs? Physiol Plant. 2005;123:153–60.

    Article  CAS  Google Scholar 

  29. Ma J, Zhang M, Xiao X, You J, Wang J, Wang T et al. Global Transcriptome Profiling of Salicornia europaea L. Shoots under NaCl Treatment. Zhou D, ed. PLoS ONE. 2013;8(6):e65877. doi:10.1371/journal.pone.0065877.

  30. Magrane M, Consortium U: UniProt Knowledgebase: a hub for protein information. Nucl Acids Res 2015;43:D204–D212.

  31. Brinker M, Brosche M, Vinocur B, Abo-Ogiala A, Fayyaz P, Janz D, et al. Linking the salt transcriptome with physiological responses of a salt-resistant Populus species as a strategy to identify genes important for stress acclimation. Plant Physiol. 2010;154:1697–709.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  32. Dang Z, Zheng L, Wang J, Gao Z, Wu S, Qi Z, et al. Transcriptomic profiling of the salt-stress response in the wild recretohalophyte Reaumuria trigyna. BMC Genomics. 2013;14:29.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  33. Fan P, Nie L, Jiang P, Feng J, Lv S, Chen X, et al. Transcriptome Analysis of Salicornia europaea under Saline Conditions Revealed the Adaptive Primary Metabolic Pathways as Early Events to Facilitate Salt Adaptation. Zhang B, ed. PLoS ONE. 2013;8(11):e80595. doi:10.1371/journal.pone.0080595.

  34. Sahu BB, Shaw BP. Isolation, identification and expression analysis of salt-induced genes in Suaeda maritima, a natural halophyte, using PCR-based suppression subtractive hybridization. BMC Plant Biology. 2009;9:69. doi:10.1186/1471-2229-9-69.

  35. Qi YC, Liu WQ, Qiu LY, Zhang SM, Ma L, Zhang H. Overexpression of glutathione S-transferase gene increases salt tolerance of arabidopsis. Russ J Plant Physiol. 2010;57:233–40.

    Article  CAS  Google Scholar 

  36. Roxas V, Smith R, Allen E, Allen R: Overexpression of glutathione S-transferase/glutathione peroxidase enhances the growth of transgenic tobacco seedlings during stress. Nat Biotechnol 1997;15:988.

  37. Hurkman WJ, Tanaka CK. Effect of salt stress on germin gene expression in barley roots. Plant Physiol. 1996;110:971–7.

    PubMed Central  CAS  PubMed  Google Scholar 

  38. Lane BG, Dunwell JM, Ray JA, Schmitt MR, Cuming AC. Germin, a protein marker of early plant development, is an oxalate oxidase. J Biol Chem. 1993;268:12239–42.

    CAS  PubMed  Google Scholar 

  39. Busov V, Johannes E, Whetten R, Sederoff R, Spiker S, Lanz-Garcia C, et al. An auxin-inducible gene from loblolly pine (Pinus taeda L.) is differentially expressed in mature and juvenile-phase shoots and encodes a putative transmembrane protein. Planta. 2004;218:916–27.

    Article  CAS  PubMed  Google Scholar 

  40. Dana MM, Pintor-Toro JA, Cubero B. Transgenic tobacco plants overexpressing chitinases of fungal origin show enhanced resistance to biotic and abiotic stress agents. Plant Physiol. 2006;142:722–30.

    Article  PubMed Central  Google Scholar 

  41. Mullen RT, Lisenbee CS, Miernyk JA, Trelease RN. Peroxisomal membrane ascorbate peroxidase is sorted to a membranous network that resembles a subdomain of the endoplasmic reticulum. The Plant Cell Online. 1999;11:2167–85.

    Article  CAS  Google Scholar 

  42. Ondrasek G. The responses of salt-affected plants to cadmium. In: Salt stress in plants. Edited by Ahmad P, Azooz MM, Prasad MNV. Springer New York; 2013: 439–463.

  43. Nedjimi B, Daoud Y. Cadmium accumulation in Atriplex halimus subsp. schweinfurthii and its influence on growth, proline, root hydraulic conductivity and nutrient uptake. Flora - Morphology, Distribution, Functional Ecology of Plants. 2009;204:316–24.

    Article  Google Scholar 

  44. Lefevre I, Marchal G, Meerts P, Correal E, Lutts S. Chloride salinity reduces cadmium accumulation by the Mediterranean halophyte species Atriplex halimus L. Environ Exp Bot. 2009;65:142–52.

    Article  CAS  Google Scholar 

  45. Guan C, Rosen ES, Boonsirichai K, Poff KL, Masson PH. The ARG1-LIKE2 gene of Arabidopsis functions in a gravity signal transduction pathway that is genetically distinct from the PGM pathway. Plant Physiol. 2003;133:100–12.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  46. Kroczyńska B, Coop NE, Miernyk JA. AtJ6, a unique J-domain protein from Arabidopsis thaliana. Plant Sci. 2000;151:19–27.

    Article  Google Scholar 

  47. Pardo J, Reddy M, Yang S, Maggio A, Huh G-H, Matsumoto T, et al. Stress signaling through Ca2+/calmodulin-dependent protein phosphatase calcineurin mediates salt adaptation in plants. Proc Natl Acad Sci U S A. 1998;95:9681–6.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  48. Galon Y, Aloni R, Nachmias D, Snir O, Feldmesser E, Scrase-Field S, et al. Calmodulin-binding transcription activator 1 mediates auxin signaling and responds to stresses in Arabidopsis. Planta. 2010;232:165–78.

    Article  CAS  PubMed  Google Scholar 

  49. Cheong YH, Kim K-N, Pandey GK, Gupta R, Grant JJ, Luan S. CBL1, a calcium sensor that differentially regulates salt, drought, and cold responses in Arabidopsis. The Plant Cell Online. 2003;15:1833–45.

    Article  CAS  Google Scholar 

  50. Nakajima S, Ito H, Tanaka R, Tanaka A. Chlorophyll b reductase plays an essential role in maturation and storability of Arabidopsis seeds. Plant Physiol. 2012;160:261–73.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  51. Lokhande V, Suprasanna P. Prospects of halophytes in understanding and managing abiotic stress tolerance. In: Ahmad P, Prasad M, editors. Environmental adaptations and stress tolerance of plants in the era of climate change. Maharashtra, India: Springer; 2012.

    Google Scholar 

  52. Martínez-Ballesta MC, Bastías E, Carvajal M. Combined effect of boron and salinity on water transport: the role of aquaporins. Plant Signal Behav. 2008;3:844–5.

    Article  Google Scholar 

  53. Jarvis DE, Ryu C, Beilstein MA, Schumaker KS. Distinct roles for SOS1 in the convergent evolution of salt tolerance in Eutrema salsugineum and Schrenkiella parvula. Mol Biol Evol. 2014;31:2094–107.

    Article  CAS  PubMed  Google Scholar 

  54. Hechenberger M, Schwappach B, Fischer WN, Frommer WB, Jentsch TJ, Steinmeyer K. A family of putative chloride channels from Arabidopsis and functional complementation of a yeast strain with a CLC gene disruption. J Biol Chem. 1996;271:33632–8.

    Article  CAS  PubMed  Google Scholar 

  55. Qiu Q, Guo Y, Dietrich M, Schumaker K, Zhu J. Regulation of SOS1, a plasma membrane Na+/H+ exchanger in Arabidopsis thaliana, by SOS2 and SOS3. Proc Natl Acad Sci U S A. 2002;99:8436–41.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  56. Hill CB, Jha D, Bacic A, Tester M, Roessner U. Characterization of ion contents and metabolic responses to salt stress of different Arabidopsis AtHKT1.1 genotypes and their parental strains. Molecular Plant. 2012;

  57. Senadheera P, Maathuis FJM. Differentially regulated kinases and phosphatases in roots may contribute to inter-cultivar difference in rice salinity tolerance. Plant Signal Behav. 2009;4:1163–5.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  58. Lee Y, Giorgi F, Lohse M, Kvederaviciute K, Klages S, Usadel B, et al. Transcriptome sequencing and microarray design for functional genomics in the extremophile Arabidopsis relative Thellungiella salsuginea (Eutrema salsugineum). BMC Genomics. 2013;14:793.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  59. Hashimoto M, Negi J, Young J, Israelsson M, Schroeder JI, Iba K. Arabidopsis HT1 kinase controls stomatal movements in response to CO2. Nat Cell Biol. 2006;8:391–7.

    Article  CAS  PubMed  Google Scholar 

  60. Ho S-L, Huang L-F, Lu C-A, He S-L, Wang C-C, Yu S-P, et al. Sugar starvation- and GA-inducible calcium-dependent protein kinase 1 feedback regulates GA biosynthesis and activates a 14-3-3 protein to confer drought tolerance in rice seedlings. Plant Mol Biol. 2013;81:347–61.

    Article  CAS  PubMed  Google Scholar 

  61. Cui MH, Yoo KS, Hyoung S, Nguyen HTK, Kim YY, Kim HJ, et al. An Arabidopsis R2R3-MYB transcription factor, AtMYB20, negatively regulates type 2C serine/threonine protein phosphatases to enhance salt tolerance. FEBS Lett. 2013;587:1773–8.

    Article  CAS  PubMed  Google Scholar 

  62. Hundertmark M, Hincha DK. LEA (Late Embryogenesis Abundant) proteins and their encoding genes in Arabidopsis thaliana. BMC Genomics. 2008;9:118.

    Article  PubMed Central  PubMed  Google Scholar 

  63. Abe H, Yamaguchi-Shinozaki K, Urao T, Iwasaki T, Hosokawa D, Shinozaki K. Role of Arabidopsis MYC and MYB homologs in drought- and abscisic acid-regulated gene expression. The Plant Cell Online. 1997;9:1859–68.

    CAS  Google Scholar 

  64. Swindell W, Huebner M, Weber A. Transcriptional profiling of Arabidopsis heat shock proteins and transcription factors reveals extensive overlap between heat and non-heat stress response pathways. BMC Genomics. 2007;8:125.

    Article  PubMed Central  PubMed  Google Scholar 

  65. Illumina. Ilumina, Inc. 2009.

  66. Robinson MD, McCarthy DJ, Smyth GK. edger: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  67. Haddad F, Baldwin K. Reverse transcription of the ribonucleic acid: the first step in RT-PCR assay. Methods Mol Biol. 2010;630:261–70.

    Article  CAS  PubMed  Google Scholar 

  68. Livak K, Schmittgen T. Analysis of relative gene expression data using real-time quantitative PCR and the 2-(delta delta C (T) method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  69. Ali Z, Zhang DY, Xu ZL, Xu L, Yi JX, He XL, et al. Uncovering the Salt Response of Soybean by Unraveling Its Wild and Cultivated Functional Genomes Using Tag Sequencing. PLoS ONE. 2012;7:e48819.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  70. Magome H, Yamaguchi S, Hanada A, Kamiya Y, Oda K: The DDF1 transcriptional activator upregulates expression of a gibberellin-deactivating gene, GA2ox7, under high-salinity stress in Arabidopsis. The Plant Journal. 2008;56:613–26.

    Article  CAS  PubMed  Google Scholar 

  71. Thomas SG, Phillips AL, Hedden P: Molecular cloning and functional expression of gibberellin 2- oxidases, multifunctional enzymes involved in gibberellin deactivation. Proceedings of the National Academy of Sciences. 1999;96:4698–703.

    Article  CAS  Google Scholar 

  72. Kore-eda S, Cushman MA, Akselrod I, Bufford D, Fredrickson M, Clark E, et al. Transcript profiling of salinity stress responses by large-scale expressed sequence tag analysis in Mesembryanthemum crystallinum. Gene. 2004;341:83–92.

    Article  PubMed  Google Scholar 

  73. Gassmann W, Hinsch ME, Staskawicz BJ: The Arabidopsis RPS4 bacterial-resistance gene is a member of the TIR-NBS-LRR family of disease-resistance genes. The Plant Journal. 1999;20:265–77.

    Article  CAS  PubMed  Google Scholar 

  74. Aguilar R, Montoya L, Sánchez de Jiménez E: Synthesis and Phosphorylation of Maize Acidic Ribosomal Proteins: Implications in Translational Regulation. Plant Physiology. 1998;116:379–85.

    Article  PubMed Central  CAS  Google Scholar 

  75. Li Q, Liu J, Tan D, Allan A, Jiang Y, Xu X, et al. A Genome-Wide Expression Profile of Salt-Responsive Genes in the Apple Rootstock Malus zumi. International Journal of Molecular Sciences. 2013;14:21053–70.

    Article  PubMed Central  PubMed  Google Scholar 

  76. Yang R, Jarvis DJ, Chen H, Beilstein M, Grimwood J, Jenkins J, et al. The reference genome of the halophytic plant Eutrema salsugineum. Frontiers in Plant Science. 2013;4:46.

    PubMed Central  CAS  PubMed  Google Scholar 

  77. Qiu Q, Ma T, Hu Q, Liu B, Wu Y, Zhou H, et al. Genome-scale transcriptome analysis of the desert poplar, Populus euphratica. Tree Physiology. 2011;31:452–61.

    Article  PubMed  Google Scholar 

  78. Droog F: Plant glutathione S-transferases, a tale of theta and tau. J Plant Growth Regul. 1997;16:95.

    Article  CAS  Google Scholar 

  79. Zhu J, Shi J, Bressan R, Hasegawa P: Expression of an Atriplex nummularia gene encoding a protein homologous to the bacterial molecular chaperone DnaJ. Plant Cell. 1993;5:341.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  80. Sleimi N, Abdelly C: Salt tolerance strategy of two fodder halophyte species: Spartina alterniflora and Suaeda fruticosa. In Cash Crop Halophytes: Recent Studies. Edited by Lieth H, Mochtchenko M: Kluwer Academic Publishers; 2003:79–86.

  81. Miller J, Arteca R, Pell E: Senescence-associated gene expression during ozone-induced leaf senescence in Arabidopsis. Plant Physiol. 1999;120:1015–23.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  82. Kudla J, Xu Q, Harter K, Gruissem W, Luan S: Genes for calcineurin B-like proteins in Arabidopsis are differentially regulated by stress signals. Proc Natl Acad Sci. 1999;96:4718.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  83. Mendoza I, Rubio F, Rodriguez-Navarro A, Pardo J: The protein phosphatase calcineurin is essential for NaCl tolerance of Saccharomyces cerevisiae. J Biol Chem. 1994;269:8792.

    CAS  PubMed  Google Scholar 

  84. Sun Y, Wang F, Wang N, Dong Y, Liu Q, Zhao L, et al. Transcriptome Exploration in Leymus chinensis under Saline-Alkaline Treatment Using 454 Pyrosequencing. PLoS One. 2013;8:e53632.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  85. Walia H, Wilson C, Zeng L, Ismail A, Condamine P, Close TJ: Genome-wide transcriptional analysis of salinity stressed japonica and indica rice genotypes during panicle initiation stage. Plant Mol Biol. 2007;63:609–23.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  86. Taji T, Sakurai T, Mochida K, Ishiwata A, Kurotani A, Totoki Y, et al. Large-scale collection and annotation of full-length enriched cDNAs from a model halophyte, Thellungiella halophila. BMC Plant Biology. 2008;8:115.

    Article  Google Scholar 

  87. Liang S, Fang L, Zhou R, Tang T, Deng S, Dong S, et al. Transcriptional Homeostasis of a Mangrove Species, Ceriops tagal, in Saline Environments, as Revealed by Microarray Analysis. PLoS One. 2012;7:e36499.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  88. Zhang J, Xie P, Lascoux M, Meagher TR, Liu J: Rapidly Evolving Genes and Stress Adaptation of Two Desert Poplars, Populus euphratica and P. pruinosa. PLoS One. 2013;8:e66370.

    Article  CAS  Google Scholar 

  89. Ouyang Y, Chen J, Xie W, Wang L, Zhang Q: Comprehensive sequence and expression profile analysis of Hsp20 gene family in rice. Plant Mol Biol. 2009;70:341–57.

    Article  CAS  PubMed  Google Scholar 

  90. Lopes-Caitar V, de Carvalho M, Darben L, Kuwahara M, Nepomuceno A, Dias W, et al. Genome-wide analysis of the Hsp20 gene family in soybean: comprehensive sequence, genomic organization and expression profile analysis under abiotic and biotic stresses. BMC Genomics. 2013;14:577.

    Article  PubMed Central  PubMed  Google Scholar 

Download references

Acknowledgements

Special thanks to Justin Page and Dr. Joshua Udall for guiding transcriptome analysis and experimental design and for Nielsen Lab members Collin Hansen, Daniel Ricks and Kevin Prier for helping with the wet lab preparations. This research has been supported by a grant from the Pakistan-U.S. Science and Technology Cooperation Program, U.S. Dept. of State and Higher Education Commission of Pakistan, and by the Department of Microbiology and Molecular Biology at Brigham Young University.

Author information

Authors and Affiliations

  1. Department of Microbiology and Molecular Biology, Brigham Young University, Provo, UT, 84602, USA

    Joann Diray-Arce & Brent L Nielsen

  2. Department of Computer Science, Brigham Young University, Provo, UT, 84602, USA

    Mark Clement

  3. Institute of Sustainable Halophyte Utilization, University of Karachi, Karachi, Pakistan

    Bilquees Gul

  4. College of Arts and Sciences, Qatar University, Doha, Qatar

    M Ajmal Khan

Authors
  1. Joann Diray-Arce
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Mark Clement
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Bilquees Gul
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. M Ajmal Khan
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Brent L Nielsen
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Brent L Nielsen.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

JDA designed and carried out experiments, including RNA isolation, cDNA library preparation, RNA-seq analysis, gene annotation, differential expression and data interpretation, and wrote the manuscript. MC provided guidance on computational analysis of the sequence data and interpretation. BG provided information and seeds for growing the plants and helped with RNA isolation and cDNA preparation. MAK provided expertise in halophyte growth and analysis, obtained funding, and assisted in preparation of the manuscript. BLN provided guidance on the project, obtained funding, and helped format and edit the manuscript. All authors read and approved the final manuscript.

Additional files

Additional file 1:

List of assembled sequences, their annotated names, gene ontology terms, associated enzymes and KEGG maps.

Additional file 2:

Summary of the number of unigenes assigned to different pathways.

Additional file 3:

Codes for bioinformatics analysis.

Additional file 4:

List of primers for qRTPCR confirmation.

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit https://creativecommons.org/licenses/by/4.0/.

The Creative Commons Public Domain Dedication waiver (https://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Diray-Arce, J., Clement, M., Gul, B. et al. Transcriptome assembly, profiling and differential gene expression analysis of the halophyte Suaeda fruticosa provides insights into salt tolerance. BMC Genomics 16, 353 (2015). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-015-1553-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-015-1553-x

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us