Abstract Salinity is one of the important abiotic stress factors that limit crop production. Common bean, Phaseolus vulgaris L., a major protein source in developing countries, is highly affected by soil salinity and the information on genes that play a role in salt tolerance is scarce. We aimed to identify differentially expressed genes (DEGs) and related pathways by comprehensive analysis of transcriptomes of both root and leaf tissues of the tolerant genotype grown under saline and control conditions in hydroponic system. We have generated a total of 158 million high-quality reads which were assembled into 83,774 all-unigenes with a mean length of 813 bp and N50 of 1,449 bp. Among the all-unigenes, 58,171 were assigned with Nr annotations after homology analyses. It was revealed that 6,422 and 4,555 all-unigenes were differentially expressed upon salt stress in leaf and root tissues respectively. Validation of the RNA-seq quantifications (RPKM values) was performed by qRT-PCR (Quantitative Reverse Transcription PCR) analyses. Enrichment analyses of DEGs based on GO and KEGG databases have shown that both leaf and root tissues regulate energy metabolism, transmembrane transport activity, and secondary metabolites to cope with salinity. A total of 2,678 putative common bean transcription factors were identified and classified under 59 transcription factor families; among them 441 were salt responsive. The data generated in this study will help in understanding the fundamentals of salt tolerance in common bean and will provide resources for functional genomic studies. Introduction Soil salinity is one of the most severe abiotic stress factors limiting the productivity of agriculture. Although most plants are glycophytes that are highly sensitive to saline environment, halophytes are plants that naturally grow under saline conditions throughout their life cycle. Salinity effects nearly 20% of all irrigated lands worldwide [29][1] and expected to reach around 50% in the near future [30][2]. A soil is considered saline if the electrical conductivity of its saturation (EC) is above 4 dS/m [31][3] which is equivalent to approximately 40 mM NaCl. As a member of grain legumes and a glycophyte crop, common bean (Phaseolus vulgaris L.) is a major source of human dietary protein, minerals, vitamins, and represents nearly half of the consumed grain legumes worldwide [32][4]. Common bean is also vital in agriculture as it forms root nodules via symbiotic associations with nitrogen fixing bacteria [33][4]. Nearly 60% (360 Mt) of the annual fresh bean production in Turkey (FAO:[34]http://faostat.fao.org/faostat) is highly dominated in Black Sea region where the soil salinity levels can reach up to 2–4 dS/m [35][5]. However, it is known that even at 1 dS/m salinity level, the productivity of common bean can be reduced up to 20% [36][6].Thus, understanding the fundamentals of salt tolerance in common bean, eventual development of improved varieties and their introduction to saline environments are imperative in agriculture. Although halophytes may use avoidance mechanisms, glycophytes tolerate salinity by minimizing ion disequilibrium and the consequent secondary effects [37][7]. In other words, tolerance mechanisms require concerted actions of mechano-receptors, ion transport channels, and secondary signal molecules to maintain ion homeostasis as well as cascades of gene activations for hormonal metabolism, signal transduction pathways, and stress responses [38][8]–[39][12]. Considering the multifactorial nature of tolerance responses, development of tolerant plants for the benefit of sustainable crop improvement still awaits accumulation of additional knowledge about the identity of components that are involved in this process. Recent developments in high-throughput approaches to study gene expression profile have emerged as an important tool to understand how plants respond to biotic and abiotic stresses. In the last few years, there have been accumulating reports on RNA-sequencing data and expression profiling on both model plants and agriculturally important crops [40][13]–[41][21] to identify genes involved in stress responses. Although, recently such high throughput transcriptome assemblies have been started in legumes [42][22]–[43][30], there are still only a handful of studies regarding the transcriptome analysis under abiotic stress conditions in these species. These studies include the effects of drought, saline-alkaline conditions and salt stress in gene expression profiling of chickpea [44][31], soybean [45][32], Medicago truncatula [46][33], and alfa alfa [47][21] respectively. However, there is no such transcriptome analysis under abiotic stress available yet for common bean. In combination with continuously generated reference genome sequences for diverse legume species [48][34]–[49][39], next generation RNA sequence analyses will provide valuable information for both identification and cloning of stress tolerance genes which can be used to improve varieties with enhanced tolerance mechanisms. In this study, we used the Illumina high-throughput RNA-sequencing platform for transcriptomic analysis of a salt tolerant common bean, Ispir genotype. We aimed to identify differentially expressed genes (DEGs) and related pathways by comprehensive analysis of data from both root and leaf tissues of the tolerant genotype grown under saline and control conditions in a hydroponic system. De novo assemblies of the sequencing data, functional annotations of unigenes, and their characterization with gene ontology and metabolic pathway analysis provided potential lists of candidate genes. Functional identification of these candidates using reverse genetics approaches in our ongoing studies will contribute to the understanding of salt tolerance mechanisms. Methods Plant growth and salt treatments The seeds of salt tolerant “Ispir” variety were kindly supplied by Prof. H. Yildiz Dasgan (Cukurova University Department of Horticulture, Adana, Turkey). The seeds were surface sterilized in a solution containing 5% (v/v) hypochlorite for 5 min and rinsed three times with distilled water. The seeds were germinated in vermiculite containing plug trays at 24/20°C cycle under a 16-h light/8-h dark photoperiod with 300 μmol m^−2 s^−1 light intensity, and 50–60% relative humidity up to the fully expanded foliage stage by daily watering with hydroponic nutrient solution containing 3 mM Ca(NO[3])[2], 1 mM MgSO[4], 0.9 mM K[2]SO[4], 0.2 mM KH[2]PO[4], 0.1 mM FeEDTA, 10 nM H[3]BO[3], 1 nM MnSO[4], 1 nM ZnSO[4], 0.1 nM CuSO[4], and 0.01 nM (NH)[6]Mo[7]O[24] [50][40]. Eight seedlings were transferred to two pots (four seedlings each) containing hydroponic nutrient solution that was replenished daily. During salt treatment, the photoperiod was also kept the same as in the germination conditions. Once plants reached trifoliate stage (five days post transfer), one pot was left as control and the other was exposed to gradually increasing NaCl concentrations. In this study, the hydroponics system was preferred to keep the nutrients and NaCl levels under strict control to achieve homogenous growth of the plants. Furthermore to minimize the risk of plasmolysis due to the osmotic shock [51][41] during salt treatments “gradual step acclimation” method was used [52][42], thus salt application was started with 50 mM in the first day, followed by 100 mM in the second day, and reached to 125 mM between days three to five in hydroponic nutrient solution. Sample collection and RNA isolation The root and the leaf tissues from salt treated and control plants were sampled at the fifth day of the salt treatment. Samples were frozen in liquid nitrogen and stored at −80°C prior to RNA extractions. Total RNAs were extracted with RNeasy Plant kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. The concentrations of RNA samples were determined by NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE). RNA quality was assessed by 1% denaturing agarose gel electrophoresis. Initially, RNA samples from root tissues of two control (RC: Root Control; RC1 and RC2) as well as two salt treated plants (RS: Root-Salt treated; RS1 and RS2) were isolated as two biological replicates, and they were pooled as described in [53]Figure 1A and 1B. Forty μg of total RNA from the two biological replicates were sequenced by Illumina HiSeqTM 2000 system (BGI, Shenzen, China). Pearson correlation coefficients between the RPKM (reads per kilobase per million reads) values of the two biological replicates were calculated as 0.99 and 0.97 for control and treated samples respectively ([54]Figure 1C and 1D). Due to the observation of high correlation within the biological replicates of root tissues, we pooled RNA samples directly from leaf tissues of the same four plants from both control (LC: Leaf Control) and treated samples as a cost effective strategy [55][43], [56][44]. Figure 1. Schematic representation of the sample pooling strategy and correlation analysis of root samples RPKM values. [57]Figure 1 [58]Open in a new tab Black lines indicate the pooling of leaf and root samples from control plants (A) and salt treated plants (B). Graphs show the correlation analysis within root control (C) and salt treated (D) samples. Nomenclature: LC: Leaf Control; LS: Leaf Salt treated; RC1 and RC2: Root Control 1 and 2; RS1 and RS2: Root Salt treated 1 and 2. cDNA library construction and Illumina sequencing RNA quality and quantity were verified using Nanodrop 1000 spectrophotometer and Agilent 2100 bioanalyzer before cDNA library generation at BGI (Shenzen, China). Total RNAs were treated with DNase I before poly (A+) mRNA enrichment with oligo dT magnetic beads. Poly (A+) RNAs were digested into 200–700 nt fragments by RNA Fragmentation Reagent, and random hexamer primed poly (A+) RNA fragments were transcribed into first-strand cDNAs. Subsequently in the presence of dNTPs, RNase H, and DNA polymerase I, the second strands were synthesized, purified using QiaQuick PCR purification kit (Qiagen, Hilden, Germany), and used for end repair single adenine nucleotide addition. The sequencing adaptors were ligated to the fragments. The paired end library constructions were finalized by size selection with agarose gel electrophoresis and selective enrichment of the cDNA fragments with PCR amplification. The libraries were sequenced on a flow cell using Illumina HiSeq2000 sequencing instrument after quality control with Agilent 2100 bioanalyzer and qPCR to detect fragment size and concentration. De novo assembly and data analyses Subsequent to adapter trimming, the raw data were filtered for reads with more than 5% ambigous bases and/or low quality reads with bases 20% of which has a Phred score less than 10. The clean reads from each library were de novo assembled into contigs with Trinity software [59][45] ([60]http://sourceforge.net/projects/trinityrnaseq/files) setting k-mer length to 25. This has generated a total of six subtranscriptomes. Contigs of the subtranscriptomes were pooled, clustered, and assembled using the Trinity software to obtain sequences that can no longer be extended on either end, and they were referred as all-unigenes. All-unigene sequences were aligned using BLASTx against NCBI non-redundant (Nr) protein, Swiss-Prot protein, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and Cluster of Orthologous Groups (COG) databases to determine sequence orientations and protein coding region predictions. Proteins with the highest ranks in BLASTx alignment results were used to predict coding region sequences. All-unigenes that could not be aligned to sequences in any of the databases mentioned above were scanned by ESTScan software [61][46] to predict sequence orientations and coding regions. For annotations, all-unigenes were searched against the Nr database using BLASTx with 10^−5 as E-value cut-off point and sequences with the highest similarities were retrieved. To obtain Gene Ontology (GO) terms regarding biological process, molecular function and cellular component [62][47] descriptions, the resulting BLASTx hits were analyzed by Blast2GO software [63][48]. The GO annotations were functionally classified by WEGO [64][49] software for gene function distributions of common bean species at macro level. BLASTx analysis against the KEGG pathway database was also performed to assign putative metabolic pathways to all-unigenes. To estimate gene expression levels, six (four from the root libraries, two from the leaf libraries) RPKM values were calculated for every all-unigene by mapping six subtranscriptomes with SOAP2 software [65][50], applying three mismatches as threshold. The clean reads mapped to more than one all-unigene were not used to calculate RPKM. Corrections for false positive and false negative errors were performed by calculating the FDR (false discovery rate) values [66][51]. The DEGs were selected by using a FDR-value ≤0.001 and the absolute ratio of log[2] (RPKM-tr/RPKM-cont) ≥1 as threshold values. The GO terms and the KEGG pathways that were enriched within the DEGs were identified by publicly available agriGO [67][52] and FatiGO [68][53] software respectively. To analyze the functional significance of enriched terms hypergeometric tests were employed by using the common bean transcriptome as background, setting the FDR and the Adjusted P-values lower than 0.05 for the agriGO and FatiGO software respectively. qRT-PCR analyses Quantitative Reverse Transcription PCR (qRT-PCR) analyses were performed for 43 unigenes ([69]Table S1) using the insulin degrading enzyme (IDE, Unigene29213) [GenBank: [70]FE702602.1] and the actin-11 (Act-11, CL442.Contig3) [GenBank: [71]CV529679.1] genes of common bean as stably expressed internal references under salt stress [72][54].