Seeds of the self-pollinated wild H. spontaneum were collected from a location in Rafah, North Sinai, Egypt (31.313559, 34.205973) near the Mediterranean Sea (∼38 g/L), where no specific permission is required. Besides, no endangered or protected plant or animal species are grown in this location. Seeds were germinated in the greenhouse in trays filled with potting mix consisting of vermiculite:perlite (1:1) and grown at 14 h of light per day, 80% humidity and 22 °C for two weeks and watered with half-strength Hoagland solution [22]. Seedlings were then salt stressed (500 mM NaCl or 29.22 g/L in half-strength Hoagland solution) at 0, 2, 12 and 24 h time intervals. Leaves of individual plants in three replicates were harvested at each time point except at time point 0 where only two replicates were gathered. All tissues were flash-frozen in liquid nitrogen and stored at −80 °C.

Flash-frozen leaf material from individual plants was crushed into a fine powder in a microcentrifuge tube using a sterilized metal rod. Total RNAs were extracted from similar-sized leaf samples collected from emergent leaves using Trizol (Invitrogen, Life Tech, Grand Island, NY, USA) and treated with RNase-free DNase (Promega Corporation, Madison, WI, USA) in the presence of 1 U/μL of RNasin^® Plus RNase Inhibitor (Promega) for 2 h at 37 °C. RNAs were quantified and 30 μg (400 ng/μL) was used for RNA-Seq. To test for the presence of DNA contamination in RNA samples, the actin gene was amplified by PCR of the original RNA samples. Purified RNA samples were shipped to Beijing Genomics Institute (BGI), Shenzhen, China in three replicates of each treatment for deep sequencing and generation of datasets (at least 100 million reads per sample).

Filtered reads were aligned with up to two mismatches to the cultivated barley genome as the reference after downloading the H. vulgare unigene transcript sequences from the NCBI database (http://www.ncbi.nlm.nih.gov/unigene). This database consisted of 26,941 transcripts including those annotated as complete and partial CDSs. RSEM v1.1.6, an RNA-Seq quantification tool, was used to estimate the relative abundances and expected read counts for the transcripts. By default, RSEM uses the Bowtie aligner (Bowtie v0.12.1) to map the reads against the transcripts. Transcript quantification of the reference-aligned reads was performed with RSEM, which allowed for the assessment of transcript abundances based on the mapping of RNA-Seq reads to the assembled transcriptome.

Expected read counts were used as input to differential expression analysis by EdgeR (version 3.0.0, R version 2.1.5). Because we had three biological replicates per time point, the median of these values was used as the common dispersion factor for differential expression (DE) analyses. The remaining unmapped sequences were re-aligned against the contigs collectively assembled de novo using the Trinity RNA-Seq Assembly package (r2013-02-25) from total unmapped sequences of all treatments and replicates. Trinity was selected for transcriptome assembly based on recent studies that showed that it performs better than other available methods [23]. DE transcripts were annotated using blast-2-GO software (version 2.3.5, http://www.blast2go.org/). Blastx was performed against the NCBI non-redundant protein database with an E-value cut off of 1 e⁻⁵. GO terms were obtained for mapped and unmapped barley transcripts with the default parameters.

To identify clusters with functional enrichment, we determined a significant Pearson correlation through permutation analysis [24]. The resulting clusters were refined by visual inspection and analyzed for GO term enrichment using Blast2GO (http://www.blast2go.org/). We also clustered the RPKM data to provide a representation of absolute abundance of the transcripts.

Six transcripts were randomly selected for validating the RNA-Seq data by real time PCR with the actin gene as the reference [25]. Primers were designed using Netprimer software (http://www.premierbiosoft.com/netprimer/index.html) with the following criteria: length ∼20 bases, GC content ∼50%, minimal secondary structures, comparable annealing temperatures (55 °C) of the primer pairs, and PCR products of ∼500 bp. Total RNAs were extracted from individual plants salt-stressed for 2, 12 and 24 h and control plants. Extraction was done in three replicates and RNAs from each treatment were then bulked. Expression levels of transcripts were detected by real time PCR using the Agilent Mx3000P qPCR Systems (Agilent technology, USA). First-strand cDNA was synthesized using 1 μg of total RNA, 0.5 μg of reverse primers of each gene (Table S1) and Superscript II reverse transcriptase (Invitrogen). All cDNA-synthesized samples were diluted (1:10) prior to amplification. The reaction (25 μL) components were 12.5 μL Maxima™ SYBR Green/ROX qPCR master mix, 0.2 μM of each gene forward and reverse primers (Table S1), and PCR-grade water was added up to 22.5 μL. Finally, 2.5 μL of diluted cDNA template were added to the reaction mix. Forty PCR cycles for each gene product included denaturation at 94 °C for 15 s, annealing at 55 °C for 30 s, and extension at 72 °C for 45 s. Amplification for each sample was carried out in triplicate along with a no-template control (NTC, PCR-grade water). Data was collected and amplification plots of ΔRn versus cycle number were generated for analysis. Calculations were made to detect the expression level of each gene under a given treatment relative to its expression under control condition.

Sequencing of cDNA samples yielded between 103 and 115 million reads corresponding to over 10 billion nucleotides per sample (Table 1). The raw sequencing reads were deposited in the Small Read Archive (SRA) at GenBank (accession number SRP032854). Between 74.8 and 80.3% of the wild barley reads could be mapped to the cultivated barley reference genome; the remaining 20–25% matched no sequences (Table 1). The percentage of transcripts of wild barley found in H. vulgare unigene database (unigene-mapped) ranged from 77.4 to 81.7%. All mapped sequences were in exonic regions of the genome. The unmapped wild barley reads of all treatments and replicates were assembled de novo and the resulting contigs were used as a new reference genome. Then, the unmapped wild barley reads were aligned to this new reference genome and results indicated that between 94.3 and 95.3% of them were de novo-mapped to the new reference. The total number of transcripts generated from alignment with both reference genomes ranged from 54,572 to 59,353. The percentage of transcripts generated from alignment with the wild barley genome contigs ranged from 60.9 to 63.9%.

RNA-Seq data was used to detect the differential expression (DE) of previously annotated barley transcripts, as well as novel transcripts uncovered in this study. To statistically obtain confirmation of the differences in gene expression across treatment time, RPKM-derived read counts were compared using a likelihood ratio test [21]. Statistical analysis was reliable when applied to genes with an RPKM value ≥ 2. To determine DE transcripts, a two-fold (or greater) change in expression and false discovery rate (FDR) of 10⁻³ or less was required. The resulting number of DE transcripts was 9277 (Fig. 1a), 3861 of which represent DE transcripts from the alignment with H. vulgare reference genome in the unigene database (Fig. 1b). The remaining DE transcripts (5416) resulted from the alignment with the generated wild barley reference genome (Fig. 1c). The number of unigene-mapped DE transcripts with no blast hits was less than 100, while that of de novo-mapped transcripts was over 2300 (Fig. 1b and c, respectively).

Annotated unigene- and de novo-mapped transcripts (5100) were utilized in generating the expression clusters across time of salt stress treatment. Expression profiles of the DE transcripts were determined by a cluster analysis based on the k-means method using Pearson's correlation distance so that the similarity in relative change for each transcript or among transcripts across time of salt treatment was determined. These data were then subjected to hierarchical clustering using the Pearson correlation as the distance metric. Fig. 2 shows the expression clusters for the DE transcripts mapped on the H. vulgare reference genome (Fig. 2a) and those mapped on the de novo-assembled wild barley reference genome using the unigene unmapped reads (Fig. 2b). Two-dimensional hierarchical clustering classified DE profiles into nine expression clusters according to the similarity of their expression profiles. Visual inspection of these expression groups suggested diverse and complex patterns of regulation. Four out of the nine expression clusters were selected for further comparisons across the time of salt tolerance treatment (Fig. 3). In general, heterogeneity and redundancy were two significant characteristics for selection. Another criterion for selection was the importance of the expression pattern. We did not have a cluster for genes that were upregulated at 2, 12 and 24 h. All upregulated genes at the 2 h time point returned to the 0 h expression level either at 12 or 24 h time points (Fig. 3a–b). We selected both expression patterns as two clusters of upregulation. The third cluster includes transcripts that were downregulated at 2, 12 and 24 h time points (Fig. 3d), and the fourth represents downregulated transcripts at 2 and 12 h time points then returned to the 0 h expression level at 24 h time point (Fig. 3c).

Real time PCR of six randomly selected DE transcripts from the mRNA-Seq data resulted in successful amplification of the bands of expected sizes (Fig. S1). Quantification of the band intensities in relation to the actin control supported the direction of change of expression as detected by mRNA-Seq for these transcripts in which expression pattern of four of them fit within cluster 2, while those of the other two fit within cluster 3.

Differentially expressed genes were assigned to functional categories using blast2GO (http://www.blast2go.org/), which provided a valuable resource for detecting specific processes, functions, and pathways during salt stress in wild barley. The results indicated that 3184, 4189 and 4106 transcripts were assigned to the three main categories; “biological process”, “cellular component” and “molecular function”, respectively (Tables 2–4). Most transcripts were assigned to the cluster 3 pattern of gene expression for downregulated transcripts across time of treatment followed by cluster 2 for upregulated transcripts at 2 and 12 h time points (Fig. 3). The number of transcripts recovered by gene ontology analysis for the three main categories was higher than the total number of DE transcripts, which likely indicates that some transcripts were assigned to more than one category. The same criterion was observed for the number of subgroups in a functional group. The numbers of functional groups for the three main categories were 21, 8 and 12 transcripts, respectively (Tables 2–4). The numbers of subgroups within groups of the three categories were 71, 12 and 49 transcripts, respectively. In the three categories, numbers of groups for cluster(s) with more than 1000 transcripts were 2, 3 and 2, respectively. These groups were “cellular process” and “metabolic process” for the “biological process” category, “membrane”, “organelle” and “cell” for the “cellular component” category, and “binding” and “catalytic activity” for the “molecular function” category. The numbers of subgroups for cluster(s) with more than 1000 transcripts were 3, 3 and 0, respectively. These subgroups were “primary metabolic process”, “macromolecule metabolic process” and “cellular metabolic process” for the “biological process category”, and “membrane-bounded organelle”, “organelle part” and “cell part” for the “cellular component” category (Tables 2–4). There were some groups and subgroups of the three main categories in which DE transcripts were either upregulated (clusters 1 and/or 2) or downregulated (clusters 3 and/or 4). The most important components were “cell–cell signaling” (cluster 2) and “cell death” (clusters 3 and 4) subgroups of the “biological process” category and the “electron carrier activity”, “superoxide dismutase activity”, “electron transporter transferring electrons within the cyclic and noncyclic electron transport pathways of photosynthesis activity”, “electron transporter transferring electrons within cytochrome b6/f complex of photosystem II activity” and “structural constituent of cytoskeleton” (clusters 1 and 2) subgroups of the “molecular function” category. There were six subgroups in the “response to stimulus” group that demonstrated both upregulation and downregulation of many different transcripts under salt stress. These subgroups were “cellular response to stimulus”, “response to abiotic stimulus”, “response to stress”, “response to external stimulus”, “response to endogenous stimulus”, and “response to chemical stimulus”.

Annotated unigene-mapped and de novo-mapped transcripts under subgroups of “electron carrier activity” and “response to abiotic stimulus” that were upregulated after 12 h of salt stress exposure with fold change (FC) of ≥4 (Table S2) fit within the expression pattern cluster 2. These subgroups were selected because they included the most important highly upregulated transcripts related to salt stress. There were thousands of GO hits whose expression was altered ≥4 fold under salt stress (Table S2). The number of GO hits under the selected “electron carrier activity” subgroup was 44, 13 of which were unigene-mapped and 31 were de novo-mapped. These GO hits represented nine genes/gene families, one unigene-mapped and eight de novo-mapped. The number of GO hits under “response to abiotic stimulus” subgroup was 866, 432 unigene-mapped, 99 de novo-mapped and 335 unigene/de novo-mapped. These GO hits represented 88 genes/gene families, 55 unigene-mapped, 15 de novo-mapped and 18 unigene/de novo-mapped.

Many cytochrome p450 proteins detected in salt-stressed wild barley are likely involved in electron transfer chains as a mechanism of salt tolerance regulation. p450 usually acts as a terminal oxidase in electron-transfer chains under salt stress with a number of fundamental redox domains, e.g., FAD-containing flavoproteins, ferredoxins [27]. Besides cytochrome p450, we detected a number of salt-regulated transcripts encoding electron transfer flavoproteins and ferredoxins. Flavoproteins function in detoxifying salts in the plant cell [28,29], while ferridoxins were recently reported for their role in abiotic stress signaling [29]. Other detected transcripts encoding vacuolar cation as well as cation proton exchanger in wild barley may act as key factors in the sequestration of sodium (Na⁺) into vacuoles to avert ion toxicity in the cytosol of plants under salinity stress. Upon influx of N+ into the cell, different ATPases (PM-ATPases, V-ATPases, and V/H-ATPases) are activated [30]. This results in Na+ efflux into the outer rhizosphere by PM Na⁺/H⁺ antiporters and/or influx into vacuoles by tonoplast Na⁺/H⁺ leading to cellular ion homeostasis, consequently salt tolerance [31]. The cytoplasmic domain of the vacuolar H+-ATPase (V-ATPase), whose transcripts were also salt-regulated in wild barley, were reported to present in a SOS (salt overly sensitive)-containing protein complex [32] with a key role in regulating ion transport under salt stress. Regulation of V-ATPase activity represents an additional key function of SOS₂ in the coordination of ion transport changes during salt stress, thus promoting salt tolerance [32]. Upregulation of transcripts encoding plasma membrane (PM)-ATPase in wild barley complements those of other ATPases, as it supported the occurrence of an electrochemical gradient generated under salt stress in the intertidal C₄ grass Spartina patens [33]. Although wild barley plants in the present study were stressed by NaCl, an additional transcript encoding calcium-transporting ATPase was upregulated. This enzyme transport protein in the plasma membrane serves to remove calcium (Ca²⁺) from the cell [34].

Flavonoid biosynthesis seems to be one of the biological processes regulating salt tolerance in wild barley. One upregulated transcript, major latex protein gene (MLP) (Table S2), was previously studied in cotton (Gossypium hirsutum) and expression in roots was induced by salt stress [35]. The Gh-MLP promoter contains potential cis-acting elements for response to salt stress and fungal elicitors. Results of RT-PCR showed that expression of Gh-MLP in Arabidopsis is rapidly induced by NaCl, and induction was maintained over 72 h [35]. In our case, expression of this gene dropped back to the control level after 24 h of salt stress exposure. Gh-MLP-transgenic A. thaliana plants showed enhanced salt stress tolerance and transgenic plants allowed seeds to germinate normally after treatment with 75 mM NaCl. Total flavonoid was also enhanced in transgenic Arabidopsis as compared to the control, suggesting that Gh-MLP might also be involved in altering flavonoid content under salt stress [35]. Two wild barley transcripts coding for chalcone synthase (CHS) and two others coding for chalcone isomerase (CHI) were reported to interact under salt stress [36]. CHS, which belongs to a family of polyketide synthase enzymes (PKS), is ubiquitous in higher plants and serves as the initial step for flavonoid biosynthesis [36]. CHI was recently reported to enhance salt tolerance in salt-sensitive yeasts [37]. In the salt-tolerant Millettia pinnata, the level of transcripts involved in flavonoid biosynthesis showed the most remarkable change under salt stress [37]. A transcript for isoflavone reductase-like enzyme was also upregulated in the present study. The enzyme is a key in the isoflavonoid phytoalexin biosynthesis pathway [38] and over expression in transgenic rice (Oryza sativa) reportedly conferred resistance to reactive oxygen species (ROS) stress [39]. Two additional wild barley transcripts encoding transcription factors of the basic helix-loop-helix (bHLH) superfamily are involved in a wide range of growth and developmental signaling pathways, including abscisic acid signaling [40], flavonoid biosynthesis [41] and abiotic stress [42].

It is well-known that cytosolic superoxide dismutase enzymes, either cytSODs (Cu,ZnSOD) or mitSODs (MnSOD), act as antioxidants and protect cellular components during salt stress from being oxidized by reactive oxygen species (ROS) [43–45]. Our results suggest that this mechanism may also be playing a role in salt stress tolerance of wild barley as three transcripts for superoxide dismutase (SOD) activity were highly expressed, while two were less expressed (FC ≤ 4) under salt stress. Other wild barley transcripts encoding electron transporter iron ion binding proteins function in iron-sulfur cluster binding. In several studies, iron ion binding proteins function in detoxification of ROS under salt stress [46–48].

A highly expressed dehydrin (DHN) transcript (FC = 8.44) was also detected. DHNs, or group 2 late embryogenesis abundant (LEA) proteins, have several roles in protecting the plant cell from dehydration stress. One of them exhibits ROS scavenging [49] mediated by the interactions between their amino acid residue and ROS species (e.g., superoxide anion radical O₂⁻; singlet oxygen ¹O₂; hydroxyl radical HO⁻; Hydrogen peroxide H₂O₂). DHNs also function as antioxidants [50], ion sequestrants [51], and metal ion transporters in plant phloem sap [52]. Under reduced hydration, the K-segments of DHNs adopt α-helical conformation [53]. The amphipathic α-helices can interact with partly dehydrated surfaces of various other proteins protecting them from further loss of the water envelope. DHNs also act as “space-fillers” in which they participate in keeping the original, non-harmful distances among different intracellular complexes that helps maintain the original cell volume, thus preventing cellular collapse [49].

Ethylene production is suggested to be a mediator of the stress responses in wild barley. A detected transcript, encoding the ethylene-forming enzyme (e.g., 1-aminocyclopropane-1-carboxylate [ACC] oxidase or ACO), was shown earlier to be indirectly induced under salt stress [54] as it relies on the expression of another set of genes, ethylene receptor genes (such as NTHK1 in tobacco, see Cao et al. [55]). Recent studies have demonstrated that salt tolerance was induced by exogenous 1-aminocyclopropane-1-carboxylate in Arabidopsis [56,57]. More recently, Li et al. [58] indicated that ethylene production and activities of ACO in cucumber (Cucumis sativus) seedlings were increased significantly under salt stress (75 mM). Another detected transcript encoding ethylene-responsive element binding factor (ERF) was considered crucial in earlier studies of cotton under stress [59]. Over expression of sugarcane and soybean ERFs in tobacco conferred tolerance to dehydration [60] and high salinity [61] (200 mM NaCl) stresses.

Transcripts encoding two serine threonine-protein kinases (STKs), SpkG and SkpC, were salt-regulated in wild barley. These STKs were previously reported to regulate various cellular functions including stress responses [62,63]. Nonetheless, Liang et al. [64] indicated that growth characteristics of a spkC mutant of the unicellular cyanobacterium Synechocystis were not affected under high salt stress conditions, while, growth of a spkG mutant was completely impaired. The spkG gene plays an essential role in sensing high salt signal directly, rather than mediating signals among other kinases [64,65]. Two transcripts encoding mitogen-activated protein kinases (MAPK) and MAP2K seem to complement the effect of SpkG in wild barley as the MAKP superfamily is part of the serine/threonine kinases. This superfamily is a key player in some of the most essential roles in plant signaling networks, and is tolerant to a variety of stresses including drought and salinity [66].

Salt stress, like other stresses, results in aggregation of cytoplasmic proteins. The present study is the first to indicate the possible role of a transcript encoding ATP-dependent Clp protease adaptor protein (ClpS) in plants with FC of 9.4. This de novo-mapped transcript was previously recovered in rice (NCBI database, http://www.uniprot.org/uniprot/Q0JNQ7) but information on its response to salt stress was not reported. Expression of ClpS in a metagenomic clone recovered [26] was shown to be involved in increased recognition of aggregated protein for refolding or degradation by the ClpA protease (ClpAP) complex resulting in salt stress tolerance. The Clp family has been shown to act as a molecular chaperone in bacteria [12]. These proteins can reverse protein aggregation resulting from salt stress and also play a role in ATP-dependent degradation of polypeptide chains under salt stress [67].

We can conclude that mRNA-Seq is an efficient high-throughput method for analyzing the wide diversity of genes expressed under salt stress. This technology is a very valuable tool to enhance our understanding of the genetics underlying salt stress tolerance mechanisms in plants. In the present study, a valuable dataset of thousands of DE transcripts was detected, some of which are involved in novel biological processes regulating salt tolerance. Some of the most likely candidates involved in salt stress in wild barley are genes involved in electron transport and exchanger mechanisms, flavonoid biosynthesis, superoxide dismutase (SOD) activity, ethylene production, signaling network, and protein refolding. The results of these comparisons can be utilized to improve salt stress tolerance in cultivated barley as well as other important cereal crops.