Transcriptome profiling of the salt-stress response in Triticum aestivum cv. Kharchia Local

Kharchia Local wheat variety is an Indian salt tolerant land race known for its tolerance to salinity. However, there is a lack of detailed information regarding molecular mechanism imparting tolerance to high salinity in this bread wheat. In the present study, differential root transcriptome analysis identifying salt stress responsive gene networks and functional annotation under salt stress in Kharchia Local was performed. A total of 453,882 reads were obtained after quality filtering, using Roche 454-GS FLX Titanium sequencing technology. From these reads 22,241 ESTs were generated out of which, 17,911 unigenes were obtained. A total of 14,898 unigenes were annotated against nr protein database. Seventy seven transcription factors families in 826 unigenes and 11,002 SSRs in 6,939 unigenes were identified. Kyoto Encyclopedia of Genes and Genomes database identified 310 metabolic pathways. The expression pattern of few selected genes was compared during the time course of salt stress treatment between salt-tolerant (Kharchia Local) and susceptible (HD2687). The transcriptome data is the first report, which offers an insight into the mechanisms and genes involved in salt tolerance. This information can be used to improve salt tolerance in elite wheat cultivars and to develop tolerant germplasm for other cereal crops.

Scientific RepoRts | 6:27752 | DOI: 10.1038/srep27752 more indepth analysis. Hence, there is a need to develop a clear understanding of salt stress tolerance mechanism of plants and identify the mechanism involved in imparting salt tolerance.
It is well known that hexaploid bread wheat shows higher salt tolerance than its tetraploid progenitor, Triticum turgidum 14 . The higher ploidy level and complex genome composition of hexaploid wheat, provides it greater physiological and ecological plasticity than its tetraploid and diploid progenitors 15 . The salt-tolerant trait in wheat cultivars grown in India, has been derived from land race Kharchia local, which is a selection from farmer's field of salinity affected Kharchi-Pali area of Rajasthan 16 . Therefore, the salt tolerant Kharchia local landrace is an excellent genetic resource for deciphering mechanism of salt tolerance in wheat (Triticum aestivum).
RNA-sequencing techniques have proven to be advantageous and economical for transcriptomic studies both in model and non-model plant species, as they neither require previously annotated genome, nor pre-synthesized nucleotide as probes. These techniques are also not limited by Expressed Sequence Tag (EST) availability 17 . The 454-FLX massively parallel DNA sequencing platform is a widely used next generation sequencing technology. It has enabled studies of the global transcriptome of even non-model plant species 18 . Large number of studies on the high throughput sequencing of plant transcriptome in many model and non model species, such us maize, grapevine, Taxus, eucalyptus, olive and cucumber etc. have been undertaken [19][20][21][22][23] . Comparative global transcriptome analysis is a powerful approach for discovering the molecular basis underlying specific biological events [24][25][26] . This approach not only allows us to map change in gene expression under specific conditions, but also to detects unique transcripts [27][28][29] .
In the present investigation, we studied the expression patterns of genes involved in salt tolerance of Kharchia Local. A comparison of the expression pattern of some of the selected genes was carried out in Kharchia Local (salt tolerant) and HD2687 (susceptible) cultivar of wheat. To our knowledge, this is the first attempt to analyse T. aestivum cv. Kharchia Local root transcriptome under salt stress. The data generated from this study, will serve as a blueprint of gene expression profile in roots of wheat under salt stress. Comparative root transcriptome profiling under salt stress, will also help to unravel the network of differentially expressed genes in the salt-tolerant wheat cultivar.

Results
Physiological and Biochemical Analysis. Physiological and biochemical analysis i.e. Relative Water content (RWC), Membrane Stability Index (MSI), Chlorophyll content (CHL), Carotenoid (CAR) and Proline content were carried out in control and salt stressed seedlings of wheat cv. Kharchia local. We observed decrease in RWC, MSI, CHL and CAR content under salt stress (250 mM NaCl) in comparison to control samples ( Supplementary Fig. S1a-d). An increase in proline and total soluble sugars content under salt stress was observed ( Supplementary Fig. S1e,f).
Sequencing and de novo transcriptome assembly of T. aestivum. In order to study gene expression pattern in response to salt stress, total RNA was extracted from root samples of T. aestivum control (CWR) and salt treated (TWR) seedlings. There were three biological replicates and three technical replicates of each, from which equimolar amount of RNA was pooled for cDNA synthesis. High throughput sequencing was performed using Roche 454 GS FLX Titanium sequencing platform. This sequencing run produced a total of 106,649 and 347,233 raw reads in CWR and TWR samples, respectively. After removal of low quality, adapter and barcode sequences a total of 367,022 (CWR-98,245; TWR-268,777) high quality reads totaling 222,600,000 bp were obtained. An overview of the sequencing and assembly process is presented in Table 1. De novo assembly of these high quality cleaned reads generated 22,531 ESTs. Using CD-HIT (V.4.6.1) software 30 , the ESTs were further assembled into 17,911 unigenes, with a minimum unigene size of 201 bp, a maximum size of 3,602 bp and an average length of 580 bp. The size distribution of the unigenes is shown in Supplementary Fig. S2a. Among these unigenes, 8,573 (47.86%) were longer than 500 bp, out of which 1,950 (22.75%) were greater than 1,000 bp. GC content gives important indication about the genes and genomic composition including evolution, structure (intron size and number) and gene regulation. It is also an indicator of stability of DNA. The average GC content of T. aestivum unigenes was 51.4% ( Supplementary Fig. S2b), which is in range of GC levels of coding sequences in monocots 31 . The transcriptome dataset generated (raw data) from the present study has been submitted at the  Table S1, Fig. 1a). Among the aligned unigenes, 54.32% had an E-value of less than 1.0 E −50 and showed very strong homology with the gene sequence in the database. The remaining 45.68% unigenes had an E-value ranging from 1.0 E −4 to 1.0 E −50 (Fig. 1b). The similarity distribution showed that 51.45% of the aligned unigenes had homology more than 90% whereas, 48.15% showed homology between 50 to 90% and only 0.40% had lower than 50% (Fig. 1c) Fig. 1d. Sequence homology based on GO classification using Blast2GO tool revealed that out of all the assembled unigenes 12,199 were summarized under three main GO categories, including 32 functional groups (Fig. 1a). A total of 71,516 GO assignments were obtained, of these, 48.66% comprised the largest category -biological processes, followed by cellular component (26.95%) and molecular functions (24.39%) at level 2. In "biological process" category, the majority of unigenes were involved in "metabolic process" (24.04%), "cellular process" (18.82%) and "single-organism process" (14.28%). With respect to "cellular component", "cell" (42.67%), "organelle" (32.44%) and "macromolecular complex" (9.83%) were the dominant groups. Under "molecular function", the top three categories were "binding" (43.51%), "catalytic activity" (42.60%) and "transporter activity" (5.60%) (Supplementary Table S2, Fig. 2).
For a better understanding of the interactions of the putative proteins and biological functions obtained from 17,911 unigenes, a single-directional BLAST search against KEGG protein database 35 was performed. This is an alternative approach to categorize genes functions with emphasis on the biochemical pathways. A total of 6,143 unigenes had KEGG orthologous number in our study, which in turn were involved in 310 different pathways (Fig. 1a). Summary of the sequences involved is described in Supplementary Table S3. The five largest pathway groups were "ribosome" (ko03010), "spliceosome" (ko3040), "carbon metabolism" (ko01200), "oxidative phosphorylation" (ko00190) and "RNA transport" (ko03013).

Identification of Transcription Factors. Transcription factors are important upstream regulatory protein
which plays a crucial role in various developmental processes and in responses to abiotic and biotic stresses. In the present study, we identified a total of 826 unigenes, representing 4.61% of the transcriptome classified into 77 putative transcription factors (TF) families (Supplementary Table S4). Among the 77 TF families, "WD40-like" represented the most abundant category comprising of 139 (16.82%) unigenes, followed by 107 (12.95%) unigenes representing "C2H2" and 39 (4.72%) unigenes representing "MYB-HB-like" transcription factor.

Quantification of T. aestivum transcripts and analysis of Differentially Expressed Genes (DEGs) in T. aestivum transcriptome.
To compare the expression level of unigenes in the CWR and TWR libraries, the number of clean reads were compared with assembled unigenes. Mapping of all the reads onto the non-redundant set of T. aestivum unigenes revealed that the number of reads corresponding to each unigene ranged from 3.2 to 323.73 for CWR (Reads Per Kilobase of Exon Per Million Fragments Mapped) and from 4.74 to 224.675 for TWR library, respectively, indicating a very wide range of expression levels of wheat unigenes. It also indicated that very low expressed unigenes were also represented in the present assembly. Among all the assembled unigenes, 4,531 unigenes were in TWR. Out of these 2,036 unigenes were unique in salt treated seedlings. Out of 4,666 unigenes in CWR, 2,171 unigenes were unique in control samples, illustrating that 2,495 unigenes were expressed both in control and salt treated plants. (Supplementary Table S5; Fig. 4a,b).

Gene ontology analysis.
To study the functions of DEGs, GO terms were extracted using Blast2GO tool and subjected to GO enrichment analysis. Annotation of DEGs revealed 1,977 unigenes belonged to 30 GO groups while the remaining 518 unigenes were not classified ( Fig. 2; Supplementary Table S2). GO functional enrichment analysis of differentially expressed genes revealed that they were involved in "iron ion binding", "oxidoreductase activity", "hydrolase activity", "ion binding", "cation binding", "cell growth" and "pyrophosphate activity". Most of the unigenes up-regulated in TWR were involved in "oxidoreductase activity", "heme binding", "iron ion binding", "anatomical structural morphogenesis", "ATPase activity", "single organism process" and "transition metal  ion binding". This suggested that these functions were more active under salt stress. In contrast, the GO groups "microtubule-based process", "chromatin binding", "motor activity" and "cytoskeleton protein binding" were over-represented in the up-regulated unigenes of the CWR ( Table 3).

Validation of DEGs by quantitative real time RT-PCR.
To further validate the results from the 454 sequencing data, eight salt induced unigenes were selected randomly (based on their expression patterns) for qRT-PCR analysis of samples that were treated with 250 mM NaCl for 2, 12 and 24 hours and control in both salt tolerant (Kharchia Local) and susceptible (HD2687) wheat genotypes (Table 4, Supplementary Table S6). In the four stages of stress, the expression of the unigenes was up-regulated in both genotypes. However, expression level was higher in Kharchia Local (Fig. 5a-h). Salt responsive protein and Na + /H + antiporter were not found to be significantly expressed in the transcriptome data analysis. To represent these unigenes under salt stress we created a heatmap of RPKM-normalized unigenes from 454 sequencing (Fig. 5i).   Alternative splicing and Differential splicing analysis. Assembled transcripts are considered to be splice variants if they aligned to the same locus and transcribed from the same strand. Transcripts with low read support were filtered out as the transcripts with low read support are likely to be assembly artifacts. In total 621 genes were alternatively spliced, resulting in 6,375 transcripts (Table 5). In total, we found 50 differentially expressing splice variants (Supplementary Table S7).

Identification of Simple Sequence Repeats (SSRs).
All the unigenes were used to mine potential SSRs to develop SSR markers. In total, 11,002 SSRs were identified by SSR Locator V.1 in 6,939 (38.74%) unigenes of T. aestivum. Out of these unigenes, 1,392 unigenes contained more than one SSR (Table 6). Hexa-nucleotide repeats (69.66%) were the most abundant SSRs in the transcriptome. The second major class was tri-nucleotide (19.77%) and the remaining repeat motifs were mononucleotide (0.16%), dinucleotide (1.24%), tetranucleotide (7.24%) and pentanucleotide (1.93%) (Fig. 6a). SSRs with two tandem repeats (67.71%) were the most common tandem repeat, followed by four, three, five and six tandem repeats (Fig. 6b) Table S8). The analysis of GC content of microsatellites showed that most of the SSRs were rich in GC content, except for mono-nucleotide repeats (Fig. 6c). Focusing on tri-and hexa-nucleotide repeats, SSRs with 50-100% GC content were thrice as prevalent as those with only 0-50% GC content. SSRs can be used as functional markers in various genetic studies in the future.

Discussion
With advances in functional genomics, proteomics and physiological studies, the molecular mechanisms of salinity stress tolerance are slowly being unraveled; however, more efforts are needed for deciphering the complexity of this stress factor in plants 37 . Next-generation sequencing technology provides a powerful tool for transcriptome analysis and the de novo assembly of transcript sequences is a rapid approach to identify expressed genes in non-model organisms 18 . The main aim of this study was to decipher the molecular mechanism of salt stress tolerance in wheat and to identify important genes and complex pathways that play a critical role in response to salt stress in T. aestivum at seedling stage. Salt stress at seedling stage is very crucial as most of the crop plants are sensitive at this developmental stage. We analyzed several physiological parameters under salt stress in wheat. A significant reduction in RWC, CHL, CAR and MSI and increase in proline and total soluble sugars content under salt stress was observed in the present study ( Supplementary Fig. S1). To decode the molecular aspect of salt stress tolerance and investigate the relationship between physiological and molecular kinetics, we generated RNA-seq libraries from roots of Kharchia local variety of T. aestivum seedlings under normal growth conditions and in response to 24 h of salt stress. To our knowledge this is the first report of root transcriptome profiling of T. aestivum cv. Kharchia Local at seedling stage under salt stress. We obtained over 453,882 raw reads leading to 17,911 predicted unigenes with N50 length of 648 bp and a total length of 10.39 Mb. Only 83.18% of these unigenes are represented in public database. Based on the annotation, some of the unigenes were closely related with the plant stress functions i.e. stress tolerance function (salt responsive protein, salt overly sensitive 1, salt tolerant like protein, salt induced protein, sodium hydrogen exchanger, sodium calcium exchanger, universal stress protein, stress protein, calcineurin B like protein, CBL-interacting protein kinase family), signal transduction (calmodulin, calcium calmodulin dependent protein kinase), energy production and conversion (ATP synthase beta subunit, ATP synthase subunit d, ATP synthase delta subunit, ATP binding protein, ATP citrate synthase, vacuolar ATP synthase subunit b, vacuolar ATPase b subunit, vacuolar type H + ATPase, vacuolar proton ATPase b subunit) and inorganic ion transport (Na + /H + antiporter, vacuolar proton-inorganic pyrophosphatase, transmembrane protein, plasma membrane H + -ATPase) (Supplementary Table S1). These putative functional unigenes identified in the present study, can provide leads for future investigations and valuable information for deciphering the putative function of novel genes. However, detailed studies are needed to understand their precise roles and functions.
Out of all the assembled unigenes, 12,199 unigenes were classified into 32 functional groups with 48.66% to biological process, 24.39% to molecular function and 26.95% to cellular component. This is indicative of the tissues undergoing extensive metabolic activities. Genes involved in some important biological processes such as response to stimulus (9.14%), single organism process (14.34%) and biological regulation (6.99%) were also  identified. This was mostly in harmony with other analysis of DEG(s) under salt stress conditions (Fig. 2). These results suggest that T. aestivum may have unique genes that regulate the response to salt stress. Logically, genes encoding these functions may be more conserved across different species and are thus easier to annotate in the database. Besides GO analysis, KEGG pathway mapping based on EC numbers for assignments was also carried out for the assembled sequence, which is an alternative approach to categorize gene functions with the emphasis on biochemical pathways. Orthology assignment and mapping of the contigs to the biological pathways were performed using KEGG automatic annotation server (KAAS). According to KEGG result, 6,143 unigenes were mapped onto a total of 310 predicted metabolic pathways (Table 2). Some stress signaling pathways represented in our study were "calcium signaling pathway", "oxidative phosphorylation", "fatty acid biosynthesis", "flavonoid biosynthesis", "ABC transporters" and "biosynthesis of other secondary metabolites". At whole transcriptome level, KEGG pathway and COG analysis are very useful techniques for prediction of potential/putative genes and their functions. The predicted molecular pathways, together with COG analysis, are useful for further investigations of genes functions.   A total of 77 transcription factor family was identified in this study, out of which, WD-40 was most abundant followed by C2H2, MYB-HB-like, CCHC (Zn), PHD and bZIP. It has been already reported that WD-40 protein is expressed under salt stress in various plants such as BnSWD1 fromin Brassica napus 38 , SRWD in Oryza sativa 39 , SiWD40 in foxtail millet 40 and WD40 repeat protein in Medicago trunculata 41 . Other TFs like C2H2 have been reported to play important role in biological processes and various abiotic stresses, involving cold, drought, oxidative and salt stress in rice 42,43 ; cold and drought stress in soybean 44 and osmotic, cold and mechanical stress in poplar 45 . This suggests that C2H2-ZFPs are likely to be associated with multiple physiological processes and stress responses in this salt tolerant land race of wheat also 46 . Over-expression of OsMYB3R-2 MYB gene has been shown to increase tolerance to cold, drought and salt stress in Arabidopsis 47 . Over-expression of OsMYB48-1 in rice improved tolerance to drought and salinity stresses 48 .
To identify significant gene expression changes associated with salt stress, differentially expressed genes (DEGs) were analyzed by comparing the CWR and TWR libraries. In total, 865 unigenes were categorized as up-regulated genes and 1,630 as down-regulated in TWR as compared to CWR (Fig. 4a). Of the 2,495 differentially expressed sequences, we predicted annotations for 2,329 unigenes against nr protein database (Fig. 4b). It is  important to identify and understand the function and role of genes of unknown function which exhibited 10 fold change in expression level under salt stress. This will improve our understanding of the salt tolerance mechanism of this tolerant bread wheat cultivar Kharchia local and allow us to employ these strategies in improving tolerance to salt stress in exotic varieties of T. aestivum. The salt stress tolerance in T. aestivum cv. Kharchia Local is capable of efficiently sequestering Na + into the root cell vacuoles. This probably restricts the Na + loading into the xylem. This process is implemented by transporting ions across tonoplast 49 . In the present study, V-ATPase (V-type proton-transporting ATPase; comp4743_ c0_seq3) unigene was present at high levels in response to salt stress. V-ATPase is essential for making an electrochemical H + -gradient across tonoplasts to boost tonoplasts for efficient ion uptake into vacuoles through the tonoplast Na + /H + antiporter (NHX) 37 . In our study, we observed higher transcript levels of V-ATPase, which provide more energy to form a stronger proton gradient to impel excessive cytoplasmic Na + into vacuoles. These findings suggest that T. aestivum has tighter control over ion compartmentalization by adjusting the driving force of ion transport which is not observed in other varieties of bread wheat.

Items Numbers
On exposure of abiotic stresses to plants, such as salinity, drought and low temperatures, plants exhibits an enhanced level of reactive oxygen species (ROS), superoxide radicals (O −2 ), hydrogen peroxide (H 2 O 2 ), and hydroxyl radicals ( − OH) that can disturb cellular homeostasis resulting in oxidative damage to cellular structures which leading to cell death 37 . The unigene associated with ROS scavenging-related gene "glutathione S-transferases" (GSTs; comp2115_c0_seq1) was also observed to be up-regulated in our study. Previously, over expression of GST has been reported to improve tolerance to abiotic stress in tobacco and soybean 50,51 . Peroxidases (PODs; comp103_c0_seq1), which are known to play an important role in H 2 O 2 scavenging, was also over expressed in our study, which is in accordance with the results of salt-treated halophyte Aeluropus littoralis 52 . These results suggests that the ROS scavenging pathway enzymes play an important role in protecting cells from oxidative damage under salt stress.
The adjustment of energy metabolism under salt stress conditions is an important strategy to impart tolerance to salt stress 37 . We identified two unigenes coding for "ATP citrate synthase" (comp4786_c0_seq6) and "atp synthase" (comp4755_c0_seq2) which are related to energy metabolism. The expressions of these two unigenes were induced under salt stress. Another unigene homologous to "cytochrome c oxidase" (comp5253_c0_seq6), which is an electrostatically coupled energy transducer that contributes to the formation of ATP in aerobic life was found to be up-regulated in the present study 37 .
"Cbl-interacting protein kinase (CIPK; comp1008_c0_seq1)" and "Calcineurin b like protein (CBLs; comp8612_c0_seq1)" constitute a complex signaling network acting in diverse plant stress responses. Earlier studies in Arabidopsis and rice have shown the importance of Ca 2+ signals in re-establishing cellular ion homeostasis, indicating the significance of the calcium sensor and signaling pathways involved in salt stress signal transduction 53 . The high salt stress tolerance of Kharchia local enabled us to identify a series of gene expression changes and snapshot of underlying salt-responsive genes. The enriched dominant GO terms that were identified during salt stress included "single organism process", "metal ion binding", "heme binding" which are in accordance with analysis of DEGs under salt stress conditions in other plans 37 . This suggests that different genes usually cooperate with each other to exercise their molecular function.
To further explore, we mapped DEGs on MapMan tool (Fig. 7). Among the TF family, zinc finger, bHLH, bZIP transcription factors were up-regulated in response to salt stress. Genes mediating ROS detoxification including peroxiredoxin (POD), thioredoxin (Trx) and glutathione peroxidase (GPX) have also shown up-regulation along with several key abiotic stress responsive genes involving sulphur assimilation (APR, ATPS and AKN). Sulphur assimilation involves activation of sulfate to adenosine 5′ -phosphosulfate (APS) by ATP sulfurylase (ATPS) which is phosphorylated by APS kinase (AKN). Subsequently, APS is reduced to sulfite by APS reductase (APR). Sulfite is reduced to sulfide and incorporated into cysteine which is a precursor of glutathione (GSH). We observed significant increase in the transcript levels of APR gene in response to salt stress in our study (Fig. 7).
Generation of splice variants is an important mechanism for regulating gene expression and for increasing transcriptome plasticity and proteome diversity in eukaryotes. Presence of splice variants plays important role in various physiological processes in plants. This includes response to different abiotic biotic and stresses. There are very few studies where alternative splicing in response to stress have been studied. Large-scale or transcriptome based studies of splice variants under salt stress conditions are still relatively less. In the present study, some of the unigenes containing splice variants were showing homology with (based on the annotation with NCBI Plant RefSeq database): drought and cold related proteins, glutathione-s-transferase, β -glucosidase, salt responsive protein, Na + /H + antiporter, zinc c3hc4 type family protein, atp-citrate synthase, wd-40 repeat, calcineurin b-like protein 2, sucrose phosphate synthase (Supplementary Table S7). These results suggest that these genes played important role in imparting tolerance to salt stress. Out of 50 differentially expressing splice variants, some unigenes showed homology with cell wall associated hydrolase, udp-glucuronic acid decarboxylase 1, lipoxygenase 2, phenylalanine ammonia-lyase, and sucrose synthase. There were more than two splice variants of these unigenes in the present transcriptome of Kharchia Local. Interestingly all of these unigenes are known to play important role in salinity stress tolerance in plants [54][55][56][57][58] . Ding et al. 59 reported that in Arabidopsis there is an increase in alternate splicing frequency under salt-stress conditions 59 . Alternate splicing seems to represent an additional level of gene regulation in response to salt stress, however, further comprehensive analysis is needed in this direction. It can reveal novel insights into the role of alternate splicing in response to various abiotic stresses in plants.
Molecular markers play an important role in plant biology. These are commonly used for the analysis of plant genomes and identification of the association between genomic variation and heritable traits. In our study a total of 11,002 SSRs were identified. These were mostly hexa-nucleotide repeats and tri-nucleotide (Fig. 6a,b). To make these SSR markers useful, we employed Primer 3 to design primer pairs for each SSR. In total, 8,181 primer pairs for 11,002 (74.35%) SSRs were designed from the microsatellites with sufficient flanking sequences Scientific RepoRts | 6:27752 | DOI: 10.1038/srep27752 (Supplementary Table S9). Further validation in relation phenotypic trait is beyond the scope of this report and will be discussed in another manuscript.
In conclusion, the comparative root transcriptome of T. aestivum cv. Kharchia Local showed multiple genes and pathways involvement to gain tolerance against salt stress. Our finding enriches the presently available genomic resources and will give an impetus to research on improving the salt tolerance in important wheat varieties. The data generated in our study can be a valuable resource and support genome analysis, besides aiding in developing expression analysis platforms, identification of molecular marker identification and initiate functional and comparative genomic studies. Our comprehensive study represents a valuable contribution towards strategies for improving the salt tolerance of crop plants. Total RNA isolation and mRNA extraction. Total RNA was isolated from both the samples (control, treated) using RaFlex ™ Total RNA Isolation Kit (GeNei ™ , Bangalore, India) as per the manufacturer's protocol.

Methods
Contaminating genomic DNA was removed by DNase I (Qiagen, Hilden, Germany) treatment. The quality of the total RNA was checked by running the samples on 1.2% formaldehyde agarose gel electrophoresis under denaturing conditions. mRNA was isolated from total RNA using Oligotex ® mRNA Midi Kit (Qiagen, Hilden, Germany) as per the manufacturer's protocol. The integrity of mRNA was checked on a 1.2% agarose gel containing formaldehyde and its quantity was measured using Nano Drop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Equal quantities of high-quality mRNA from biological replicates were pooled together and were used for transcriptome study. cDNA Library construction and 454 sequencing. cDNA was synthesized from the pooled RNA samples of CWR and TWR using cDNA Synthesis System (Roche, Basel, Switzerland) as per the manufacturer's protocol using hexamer primers. The cDNA samples were sheared by nebulization to produce random fragments of approximately 300-800 bp in length. The nebulized cDNAs were recovered using QIAquick ® PCR Purification Kit (Qiagen, Hilden, Germany). Rapid Library Molecular Identifiers (RLMIDs) were ligated to the fragmented purified cDNA samples according to the standard protocol. Each MID contains a barcode sequence that is used to discriminate samples from different libraries. Library quantification was done using TBS 380 Fluorometer. The Relative Fluorescence unit of each dilution was recorded and checked against RL standard curve and the sample concentration was calculated. Agilent Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) High Sensitivity DNA chip was used to check quality of the libraries.
The quality passed samples were pooled and amplified by emulsion PCR for sequencing. The sequencing of the libraries was performed on 454 -GS FLX Titanium sequencer (454 Life Sciences, Roche, Basel, Switzerland).
Sequence data analysis and assembly. The raw reads were first pre-processed by eliminating adaptor sequences using customized Perl scripts. All the sequences smaller than 60 bases were eliminated based on the assumption that small reads might represent sequencing artifacts. At the same time Q20, Q30, GC content and sequence duplication level of clean data were calculated. Finally, high quality reads were assembled into unique putative transcripts, termed as unigenes (including contigs and singletons), using Trinity ESTs assembler. The assembly was performed using the default parameters.
Functional annotation, GO classification and metabolic pathway analysis of unigenes. Blast homology searches and sequence annotations were carried out using Blast2GO tool v.3.1.3 (http://www.blast2go. org) [32][33][34] . The assembled sequences were compared against the NCBI non-redundant (nr) protein database via Blast X, using an E-value cut-off of 1.0 E −3 . For gene ontology mapping, the program extracts the GO terms associated with homologies (GO; http://www.geneontology.org). These GO terms are assigned to query sequences, producing a broad overview of groups of genes cataloged in the transcriptome for each of three GO categories. The data presented herein represent a GO analysis at level 2, illustrating general functional categories. GO enrichment analysis was also done. This enrichment analysis was performed to evaluate the enrichment of various GO categories for the unigenes having 2 fold or above expression values in both the samples. The unigene sequences were also aligned to the COG database to predict and classify proteins. KEGG pathways were assigned to the assembled sequences using the online KEGG Automatic Annotation Server (KAAS) 60 , http://www.genome.jp/ kegg/kaas. The single-directional best hit (SBH) method was used to obtain KEGG analysis. For identification of transcription factor in Kharchia Local transcriptome, all the assembled unigenes were analyzed against the PlantTFcat online tool (http://plantgrn.noble.org/PlantTFcat). GC content of the sequences was measured using Emboss GeeCee tool (http://emboss.bioinformatics.nl/cgi-bin/emboss/geecee). Expression analysis. The RPKM method was used to calculate the unigene expression in this study 61 . The clean reads of each library were mapped to the sequences of each unigene. The significance of difference in gene expression between the control and salt treated plants was determined by using DEGseq, an R package 62 . False discovery rate (FDR) was applied to identify the threshold of the P-value in multiple tests 63 . When FDR was less than 0.05 and log ratio greater than 1 (two-fold change) between the samples the unigenes were considered as differentially expressed.
Alternative splicing Landscape and Differential splicing analysis. To analyze the various types of Alternative splicing (AS) raw reads of both the samples (CWR and TWR) were mapped to wheat reference genome released by IWGSC (http://plants.ensembl.org/Triticum_aestivum), using TOPHAT software using default parameters, which utilizes computational methods to detect variants and splicing isoforms in short reads through merging and filtering position lists from a genomic index. For differential splicing analysis, we used Cuffdiff, which categorizes individual transcripts based on their transcription start site, thus grouping together all the transcripts with a common pre-mRNA molecule 64,65 . Cuffdiff then estimates significant differences in expression of each transcript in a group relative to one another using the Jensen-Shannon divergence metric across the different samples. Differential expression calls were made using the DESeq package.
Physiological analysis. Relative water content, chlorophyll and carotenoid content, membrane stability index, proline content and total soluble sugars were measured in control and salt treated plants at seedling stage as mentioned by Sairam et al., (2012) 66 . All the experiments were performed with three biological and three technical replicates.
Validation of genes using qRT-PCR. Seeds of Kharchia Local (salt tolerant) and HD2687 (salt susceptible) seeds were grown in hydroponics as mentioned above. The 10 days old seedlings were subjected to 250 mM NaCl stress. The root samples of salt treated samples at different time points (2, 12 and 24 h) and control were taken, and total RNA was isolated using the same method mentioned above. The salt-induced expression of unigenes with potential roles in the salt stress response was chosen for validation using qRT-PCR. Three independent biological samples of each were used in the analysis. cDNA was synthesized from 2 μ g of total RNA using Superscript ® III First strand cDNA Synthesis system (Invitrogen, USA) as per the manufacturer's protocol. Actin is reported to exhibit most stable expression in different developmental stages and stresses, was used as reference for gene expression normalization. The gene specific primers are listed in Supplementary Table S8. qRT-PCR was performed using an Agilent Mx3000P ™ PCR platform and KAPA SYBR ® FAST qPCR kit Master Mix (2X) Universal (Kapa Biosystems, Woburn, MA) as per manufacturer's instructions. The relative expression levels of the selected unigenes normalized to the expression level of actin were calculated from cycle threshold values using the 2 −ΔΔCt method. This experiment was carried out in three independent biological replicate and three technical replicate of each biological replicate.