Identification and characterization of long non-coding RNAs involved in osmotic and salt stress in Medicago truncatula using genome-wide high-throughput sequencing

Long non-coding RNAs (lncRNAs) have been shown to play crucially regulatory roles in diverse biological processes involving complex mechanisms. However, information regarding the number, sequences, characteristics and potential functions of lncRNAs in plants is so far overly limited. Using high-throughput sequencing and bioinformatics analysis, we identified a total of 23,324 putative lncRNAs from control, osmotic stress- and salt stress-treated leaf and root samples of Medicago truncatula, a model legume species. Out of these lncRNAs, 7,863 and 5,561 lncRNAs were identified from osmotic stress-treated leaf and root samples, respectively. While, 7,361 and 7,874 lncRNAs were identified from salt stress-treated leaf and root samples, respectively. To reveal their potential functions, we analyzed Gene Ontology (GO) terms of genes that overlap with or are neighbors of the stress-responsive lncRNAs. Enrichments in GO terms in biological processes such as signal transduction, energy synthesis, molecule metabolism, detoxification, transcription and translation were found. LncRNAs are likely involved in regulating plant’s responses and adaptation to osmotic and salt stresses in complex regulatory networks with protein-coding genes. These findings are of importance for our understanding of the potential roles of lncRNAs in responses of plants in general and M. truncatula in particular to abiotic stresses.


Background
Non-coding RNAs (ncRNAs) are a set of RNAs that have no capacity to code for proteins. They are used to be considered as inconsequential transcriptional "noises", because of limited information for their functions [1,2]. However, this situation is being changed. Recent studies have shown that ncRNAs play important regulatory roles in numerous biological processes [3,4].
NcRNAs are grouped into small RNAs, such as micro-RNAs (miRNAs) and small interfering RNAs (siRNAs), and long non-coding RNAs (lncRNAs) according to the length [5]. LncRNAs are defined as a group of ncRNAs that have a length of more than 200 nucleotides [6]. They are usually expressed at low levels and lacking sequence similarities among species, exhibit tissue and cell-specific expression patterns, and transcripts are localized to subcellular compartments [4,7]. LncRNAs can be further grouped into sense, antisense, bidirectional, intronic and intergenic lncRNAs according to their relative locations with protein-coding genes [8]. In Arabidopsis thaliana, >30 % of lncRNAs are intergenic, and antisense lncRNAs are also abundant [9,10].
It has been shown that some lncRNAs regulate the expression of genes in a close proximity (cis-acting) or in a distance (trans-acting) in the genome via a number of mechanisms, including modifying promoter activities by nucleosome repositioning, histone modifications, DNA methylation, activating/gathering/transporting of accessory proteins, epigenetic silencing and repression [8,11,12]. Increasing evidence supports that lncRNAs play a crucial role in disease occurrence, genomic imprinting and developmental regulation in mammals [13][14][15].
In contrast to extensive studies of lncRNAs in mammals [13,14,16,17], only a few studies have been reported of the function of lncRNAs in plants [18,19]. For example, COOLAIR and COLDAIR have been identified to be associated with FLOWERING LOCUS C (FLC) in Arabidopsis. COLDAIR includes two antisense lncRNAs transcribed from the antisense strand of FLC, while COLDAIR is an intronic lncRNA transcribed from the first intron of FLC. They have been implicated in silencing and epigenetic repression of FLC to regulate flowering time [20,21]. AtIPS1 and At4 have been shown to act as target mimics of miR399 by binding and sequestering miR399 and reduce miR399-mediated cleavage of PHO2 which is important for phosphate uptake [22,23]. Genome-wide identification of lncRNAs in A. thaliana has been reported in several studies [24][25][26][27]. In rice, LDMAR has been shown to regulate photoperiod-sensitive male sterility [28]. Bioinformatics analyses reveal that 60 % of lncRNAs are precursors of small RNAs and 50 % of lncRNAs are expressed in a tissue-specific manner [29][30][31].
Medicago truncatula is a model legume widely used in genomics, genetics and physiological studies of legumes due to its small genome size and relative ease in genetic transformation [32,33]. Legumes account for one third of primary crop production in the word and are important sources of dietary proteins for human and animals [34]. In M. truncatula, Enod40 and Mt4 involved in nodulation and phosphate uptake, respectively, have been identified as lncRNAs [35,36]. Although a recent in silico analysis of lncRNAs has been conducted in M. truncatula, only limited information is presented, because only lncRNAs with poly(A) tails have been analyzed, using less finished genome sequences available at the time [37]. As most lncRNAs have no poly(A) tails and are lowly and specifically expressed [4,16], to identify a comprehensive set of lncRNAs including nonpoly(A)-tailed lncRNAs in M. truncatula, we conducted genome-wide high-throughput sequencing of six libraries prepared using complementary sequences of synthetic adaptors. Similar to other plant species, legumes are also frequently encountered adverse environments such as osmotic and salt stresses. Previous studies of molecular mechanisms underlying plant's tolerance to abiotic stresses are mainly focused on functional studies of protein-coding genes, while few studies have systemically investigated the roles of lncRNAs in osmotic and salt stress responses of plants. In the present study, we identified a comprehensive set of lncRNAs that are responsive to osmotic and salt stresses in leaves and roots of M. truncatula using high throughput sequencing of six cDNA libraries.

Physiological response to osmotic and salt stress
Materials used to construct cDNA libraries were treated by osmotic or salt stress for 5 h. Foliar osmolality was increased from 350 mOsmol kg −1 to 450 and 390 mOsmol kg −1 , after the treatments with osmotic and salt stress, respectively (Table 1). There was a significant increase in foliar Na + concentration after 5-h salt treatment (Table 1). No effects of osmotic and salt stress on concentrations of proline (Pro) and soluble sugars were detected (Table 1). These results suggest that plants under our treatment regime are at the early stage of stress-response to activate genes and their regulatory networks.

High-throughput sequencing
Six cDNA libraries were constructed using mRNA isolated from leaves and roots of M. truncatula seedlings treated with osmotic stress (OS), salt stress (SS), and control (CK) and complementary sequences of synthetic adaptors. They were sequenced by an Illumina-Solexa sequencer. The high-throughput sequencing led to more than 90,000,000 raw sequence reads. To assess the quality of RNA-seq data, each base in the reads was assigned a quality score (Q) by a phred-like algorithm using the FastQC [38]. The analysis revealed that the data are highly credible with a mean Q-value of 36 (Additional file 1: Figure S1). Of the raw reads, more than 99 % were clean reads after initial processing ( Table 2). We performed 100 bp paired-end sequencing, and led to 56.7 G raw bases and 56.6 G clean bases in total. Data are the means ± SE (n = 3). Data with "*" or "**" indicate significant different (P < 0.05 or P < 0.01) between treatments and control

Identification and characterization of lncRNAs
The clean reads were mapped to the M. truncatula genome (Mt4.0) using the TopHat [39]. Transcripts were then assembled and annotated using the Cufflinks package [40].  Table 2). In total, 23,324 unique lncRNAs were obtained in the present study (Additional file 2: Table S1). And this number was similar to that of lncRNAs in Arabidopsis and maize [30,41]. We found that these lncRNAs were more evenly distributed across the 8 chromosomes in M. truncatula with no obvious preferences of locations (Fig. 1a). According to the locations of lncRNAs in the genome, 10,426 intronic, 5,794 intergenic, 3,558 sense and 3,546 antisense lncRNAs were identified ( Fig. 1b and e). In terms of the lncRNAs' length, the majority of lncRNAs was relatively short. For example, 84.1 % of them were shorter than 1,000 nt (Fig. 1c). Interestingly, lncRNAs and mRNAs were much more abundant in roots than in leaves, given that similar amounts of raw reads were obtained for both leaf and root samples. In all libraries, more lncRNAs were detected in roots than in leaves ( Table 2). For example, 18,275 lncRNAs were identified in roots, while there were 11,501 lncRNAs in leaves under control condition (Fig. 2a). Furthermore, we found that the accumulative frequency of lncRNAs differed in leaves from that in roots. The proportion of lncRNAs with a high level of expression was more than mRNAs in leaves, but this expression pattern was in contrary in roots under the control conditions (Fig. 1d). Moreover, these patterns of expression were not altered by treatments with osmotic and salt stress (Additional file 1: Figure S2). The lack of chloroplast-derived RNAs in roots might be a possible reason for the difference between leaves and roots.
All putative lncRNAs in M. truncatula were aligned with lncRNAs in A. thaliana from NONCODE database [42]. We can only detect 140 lncRNAs that were comparable to those lncRNAs in A. thaliana, suggesting that lncRNAs are weakly conserved between the two species (Additional file 2: Table S1). Moreover, lncRNAs which were from transposons or which encoded microRNAs were marked (Additional file 2: Table S1).

Responses of lncRNAs to osmotic and salt stresses
To identify osmotic stress-and salt stress-responsive lncRNAs, the normalized expression (fragments per kilobase of exon per million fragments mapped, FPKM) of lncRNAs was compared amongst the six libraries.
LncRNAs that were responsive to osmotic and salt stresses in leaves and roots were identified by determining the P-value and false discovery rate. To verify the results from the RNA-seq experiments, 12 lncRNAs were selected to verify their expression by quantitative realtime PCR (qRT-PCR) ( Fig. 3 and Additional file 1: Figure S3). These results indicate that our transcriptomic analysis is highly reproducible and reliable, and that lncRNAs identified from the high throughput sequencing represent real transcripts.
Transcript levels of 7,863 lncRNAs in leaves and 5,561 lncRNAs in roots were detected to be changed by the osmotic stress, and 7,361 lncRNAs in leaves and 7,874 lncRNAs in roots were identified to be responsive to the salt stress. Venn diagrams showed common and specific lncRNAs, whose expression was altered in roots and leaves by osmotic and salt stresses (Fig. 2b). Some lncRNAs in leaves and roots showed different responses to osmotic and salt stresses. There were 1,783 and 2,148 lncRNAs, whose expression was changed in both leaves and roots by osmotic and salt stresses, respectively. In leaves, more than half of stress-responsive lncRNAs were common between osmotic stress (59.6 %) and salt stress (63.7 %). However, these values were decreased to 47.0 % and 33.2 % in roots, respectively. The expression levels of 471 lncRNAs were found to be changed in the four treated samples (Fig. 2b). Among the lncRNAs, whose expression was changed in responses to osmotic and salt stresses, we further classified them to upregulated and down-regulated classes (Additional file 1:    Twenty-eight and 213 lncRNAs were found to be upregulated and down-regulated, respectively, in both roots and leaves treated with osmotic and salt stresses.

Functional analysis of stress-responsive lncRNAs
Previous studies showed that lncRNAs are preferentially located in a close proximity to genes that they regulate [13,[43][44][45]. To reveal potential functions of the identified lncRNAs, we analyzed Gene Ontology (GO) terms of genes that were co-expressed and spaced by less than 100 kb with the stress-responsive lncRNAs. We detected significant enrichments (P < 0.05) of 26  . In roots, GO term enrichments were greater than those in leaves (i.e., 52 vs 37), suggesting that roots are more sensitive to osmotic and salt stresses than leaves (Additional file 1: Figure S5, Tables S4 and S5). These findings suggest that the stress-responsive lncRNAs may regulate genes involved in many biological processes, including signal transduction, energy synthesis, molecule metabolism, detoxification, transcription and translation in response to osmotic and salt stresses.
One lncRNA may regulate multiple other lncRNAs and protein-coding genes, and vice versa [4]. To unravel the relationship among lncRNAs and protein-coding RNAs which were co-expressed and spaced by less than 100 kb, putative interactive networks were constructed using Cytoscape (Fig. 5 and Additional file 1: Figure S6). About half of them had less than or equal to three nodes like networks in Fig. 5c. More complex interactive networks were also observed. For example, thirteen proteincoding genes involved in oxidation/reduction reaction, transcription, energy synthesis and signal transduction were found to be regulated by three lncRNAs in the situation of salt stress in leaves (Fig. 5a). Two transcription factors of MYB and zinc finger families were found in the network of Fig. 5b, which may activate stress-responsive genes in the downstream under osmotic stress in roots. The expression of lncRNAs in Fig. 5c has been validated in Fig. 3. TCONS_00046739 was identified as regulator of cytochrome P450 in roots under salt stress. The targets of TCONS_00100258 and TCONS_00118328 may be two transmembrane proteins in leaves under salt stress. These networks among lncRNAs and protein-coding genes may play important roles in sensing and responding to osmotic and salt stresses. The construction of putative network based on gene expression and vicinity of the lncRNAs and protein-coding genes may not be very robust due to the Under stresses, many GO terms were enriched, such as carbonate dehydratase activity (GO:0004089) and carbon utilization (GO:0015976) that are highly significant (because of the lowest P value) in leaves under osmotic and salt stresses (Fig. 4). The carbonic anhydrase gene Medtr6g006990, belonging to these two GO terms was down-regulated by these two abiotic stresses. This gene is predicted to be regulated by the lncRNA TCONS_ 00097188 located in the upstream of the coding region of Medtr6g006990 (Fig. 6a). Carbonic anhydrase catalyzing the reversible hydration of CO 2 into bicarbonate plays an essential role in the accumulation of CO 2 in the active site of rubisco [46]. Our results suggest that TCONS_00097188 may regulate photosynthesis under the abiotic stresses by regulating the expression of Medtr6g006990.
Under conditions of abiotic stresses, signal transduction networks are mobilized to cope with the stressed environment. The pathway of phospholipids metabolism has been proposed to be an important in response to a number of abiotic stresses [47]. For example, drought and salt stresses up-regulate the expression of genes encoding phosphatidylinositol-specific phospholipase C (PI-PLC), which hydrolyzes phosphatidylinositol 4,5bisphosphate to the secondary messenger molecules inositol 1,4,5-trisphosphate and diacylglycerol [47]. In the present study, the expression of a PI-PLC gene (Medtr3g069280), which belongs to GO:0004435 (Phosphatidylinositol phospholipase C activity) and GO:0007165 (Signal transduction) was up-regulated in response to osmotic and salt stresses, and the lncRNA TCONS_00047650 was expressed from the regulatory region of Medtr3g 069280 (Fig. 6b). These results suggested that TCONS_ 00047650 may regulate the expression of Medtr3g069280.
Plants under osmotic and salt stresses often display oxidative stress symptoms as indicated by marked accumulation of reactive oxygen species (ROS), which damages membrane systems. To cope with the excessive accumulation of ROS, plants mobilize antioxidant enzymes to scavenge ROS [48]. We found that the expression of Medtr7g094600 coding for glutathione peroxidase (POD) was up-regulated in roots. We identified the lncRNA TCONS_00116877 located approximately 3.9 kb upstream of the coding sequence of Medtr7g094600 (Fig. 6c). These results suggest that TCONS_00116877 may be involved in regulating plant's tolerance to the oxidative stress by modulating the expression of POD.
Effect of salinity on plant growth can be divided into ionic toxicity and osmotic stress [49]. Plants often exhibit similar tolerance mechanisms, such as altered energy synthesis, phospholipids signal transduction and detoxification to osmotic and salt stresses [47]. In addition, we found that the expression of the Na + /H + exchanger (NHX) gene Medtr1g081900 was up-regulated by the salt stress in roots. This gene codes for a vacuolar Na + /H + antiporter mediating Na + influx into the vacuoles [50]. This gene is predicted to be regulated by the lncRNA TCONS_00020253 located in the upstream of the coding region of Medtr1g081900 (Fig. 6d). These results suggest that TCONS_00020253 is likely a regulator of Medtr1g081900.

Discussion
Less than 2 % of the human genome sequences codes for proteins [51]. However, transcription is not limited to protein-coding regions [17,52]. In fact, more than 90 % of the human genome sequences are likely transcribed [17]. These non-coding transcribed sequences are from  introns, intergenic regions or the antisense strand of protein-coding genes [16]. An increasing number of studies have shown that ncRNAs play important roles in many vital biological processes, highlighting that ncRNAs are not transcriptional "noises" [4].
Studies on lncRNAs are less extensive in plants than in mammals, and those studies are mainly conducted in A. thaliana [25,26]. In addition to cereals, legumes are the most important sources for human foods and animal feeds worldwide. Moreover, legumes are unique among cultivated plants for their ability to directly utilize atmospheric nitrogen through symbiotic interactions with the soil bacteria rhizobia [32]. According to the genome sequences of M. truncatula, only about 17 % of the sequences code for proteins [33]. Previous studies in M. truncatula have been concentrated upon protein-coding sequences associated with nodulation, abiotic stresses and developmental processes [53][54][55][56]. Several recent studies have investigated functions of small RNAs involved in nodulation and abiotic stresses [57][58][59]. In this report, we show that lncRNAs are distributed in almost the entire genome of M. truncatula, suggesting that lncRNA-coding regions are much more widespread than protein-coding regions ( Fig. 1a and b). Whole genome sequencing and annotation facilitate functional studies of protein-coding genes [32,33]. Identification and characterization of the large number of lncRNAs in M. truncatula in the present study provide valuable information for functional characterization of lncRNAs in plants in general and in legumes in particular.
In the present study, the reverse transcription was made by using complementary sequences of artificial adaptors to enrich lncRNAs with or without poly(A) tails. To distinguish sense from antisense lncRNAs, strand-specific libraries were constructed and paired-end sequencing was carried out in the present study. As a result, our results can be used to identify different types of lncRNAs to facilitate functional studies. Moreover, the abundant original data (56.7 G) generated in the present study allow us to detect lncRNAs that have low expression levels. Given that the expression of lncRNAs is highly tissue-specific [30], lncRNAs from both leaves and roots of M. truncatula were sequenced and their expression patterns were compared. In addition, we also sequenced and compared the expression of proteincoding genes in both leaves and roots under control and stressed conditions. This information is useful for predicting putative targets of lncRNAs. Furthermore, we identified common and specific lncRNAs from leaves and roots treated with osmotic or salt stresses to study potential functions of lncRNAs in plant's responses to abiotic stresses. To our best knowledge, this is the first report of a comprehensive set of lncRNAs isolated from osmotic-and salt-stress treated leaf and root samples of higher plants using high-throughput sequencing. Unlike previous studies where osmotic and salt stress-responsive lincRNAs (intergenic lncRNAs) were detected in Arabidopsis [26] and Populus [60], the present study identified all types of lncRNAs involved in osmotic and salt stresses in M. truncatula by the strand-specific sequencing. Moreover, to make sure that the putative lncRNAs in this study conform to the criteria of length and proteincoding ability, the putative lncRNAs were selected to have >200 bp in length and less than -1 for the coding potential score. These strict criteria and improved methods made the identified lncRNAs with high sensitivity and selectivity.
To minimize the adverse effects of abiotic stresses, plants have evolved a suite of responsive mechanisms [49]. There are many protein-coding genes which are identified to play regulatory roles under varying abiotic stresses, such as DREB1A/CBF3, SOS1 and so on [61][62][63][64]. However, little is known of biological functions of lncRNAs in abiotic stress responses in plants. Moreover, lncRNAs are putative potent tools for plant improvement to enhance their resistance to abiotic stresses [65]. Therefore, identification of abiotic stress-responsive lncRNAs, characterization of their functions and dissection of their regulatory networks can enhance our mechanistic understanding of plant response and adaptation to stressed environment. Several recent studies have identified lncRNAs involved in biotic/abiotic stresses in plants. Fusarium oxysporum, a soil-borne plant fungal pathogen, causes the vascular wilt disease through roots in several plant species [66]. LncRNAs that are responsive to F. oxysporum have been identified by RNA-seq, and functional characterization of these lncRNAs reveals that lncRNAs are important components of the antifungal networks in A. thaliana [66]. For abiotic stress responses, 76 lncRNAs have been identified from a fulllength cDNA library of A. thaliana [25]. Of these, 22 lncRNAs have been shown involved in abiotic stress responses; overexpression of two identified lncRNAs renders plants more tolerance to salinity. However, because the full-length cDNA library was made from mRNAs with poly(A) tails, lncRNAs without poly(A) tails have not been identified in that study. In our present study, reverse transcription was made by complementary sequences of artificial adaptors, thus, lncRNAs with or without poly(A) tails were obtained. Liu et al. [26] identified 6,484 lincRNAs, of which 1,832 lincRNAs are responsive to drought, cold, salinity and/or abscisic acid. In a recent study, a total of 504 drought-responsive lincR-NAs has been detected in Populus [60]. However, in these studies, only lincRNAs, rather than all types of lncRNAs, are analyzed. In our study, all types of lncRNAs, including those of sense, antisense, intronic and intergenic lncRNAs were identified using the advanced sequencing technology and analytic methods such as strand-specific sequencing and Cuffcompare analysis.

Conclusions
In this study, we identified 23,324 putative lncRNAs from six RNA-seq libraries of M. truncatula by highthroughput sequencing, of which 11,641 and 13,087 lncRNAs are found to be responsive to osmotic stress and salt stress, respectively. Of these, 5,634 lncRNAs are found to be responsive to both osmotic and salt stress. We analyzed GO terms of genes that either overlap with or are immediate neighbors of the stress-responsive lncRNAs. We found enrichments of GO terms in many biological processes, including signal transduction, energy synthesis, molecule metabolism, detoxification, transcription and translation. Moreover, a number of complex interaction networks were constructed based on co-expression and genomic co-location of lncRNAs and protein-coding genes. These results suggest that lncRNAs are likely involved in regulating plant's responses and adaptation to osmotic and salt stresses in complex regulatory networks with protein-coding genes. These findings provide valuable information for further functional characterization of lncRNAs in responses of plants in general and M. truncatula in particular to abiotic stresses.

Plant materials and stress treatments
Seeds of Medicago truncatula ecotypes Jemalong A17 were treated with concentrated sulfuric acid for 8 min, and then thoroughly rinsed with water. After chilled at 4°C for 2 days, seeds were sown on 0.8 % agar to germinate at 25°C till the radicals being about 2 cm. The seeds were planted in the plastic buckets filled with aerated nutrient solution under controlled conditions (26°C day/22°C night, and 14-h photoperiod). The composition of fullstrength nutrient solution is: 2.5 mM KNO 3  Three-week-old seedlings were transferred into nutrient solutions containing either 265 mM mannitol or 150 mM NaCl, which had identical osmolality, for 5 h. Leaves and roots from at least ten individual plants were collected and frozen immediately in liquid nitrogen until use. At the same time, M. truncatula seedlings grown in the fullstrength solution without mannitol or NaCl were harvested and were used as control. The regimes of treatment used in this study were chosen based on previous studies [25,67].

Construction of cDNA libraries and high-throughput sequencing
To construct libraries, total RNA was extracted from leaves and roots of seedlings grown in different solutions (osmotic stress, salt stress and control) using the Trizol (Invitrogen) according to the manufacturer's protocols. Ribosome RNA of six RNA samples was removed using Ribo-Zero™ Magnetic Kit (Epicentre). Thereafter the strand-specific sequencing libraries were constructed following a previously described protocol [68]. The paired-end sequencing (2 × 100 bp) was performed on an Illumina Hiseq2000 sequencer at the LC Biotech, Hangzhou, China.

Reads mapping and transcriptome assembling
The resulting directional 100 bp paired-end reads were quality-checked with FastQC (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/), and adapter contaminations and low quality tags in the raw data were removed. Ribosome RNA data were also removed from the remaining data by alignment. Then, the clean reads from six-cDNA libraries were merged and mapped to the M. truncatula genome sequence (Mt4.0) using the spliced read aligner TopHat [39]. To construct transcriptome, the mapped reads were assembled de novo using Cufflinks [40]. All transcripts were required to be >200 bp in length.

Identification of lncRNAs
The assembled transcripts were annotated using the Cuffcompare program from the Cufflinks package [40]. According to the annotation of M. truncatula genome sequence (Mt4.0), the known protein-coding transcripts were identified. The remaining unknown transcripts were used to screen for putative lncRNAs. The transcripts smaller than 200 bp were firstly excluded. Then, the coding potential for the remaining transcripts was calculated by the Coding Potential Calculator based on quality, completeness, and sequence similarity of the open reading frame to the proteins in the protein databases [69]. A transcript was deemed to be noncoding if the coding potentials are scored to be less than −1, which suggest that this transcript has no capacity of coding for proteins.

Analysis of differential expression patterns
Expression levels of all transcripts, including putative lncRNAs and mRNAs, were quantified as FPKM using the Cuffdiff program from the Cufflinks package [40]. Differential gene expression was determined using DESeq with a P-value < 0.05 and a false discovery rate threshold of 5 % [70].

Quantitative real-time PCR (qRT-PCR)
Total RNA was isolated using RNAiso Plus reagent (TaKaRa) and treated with RNase-free DNase I (Promega). About 0.5 μg RNA was reverse-transcribed into firststrand cDNA with PrimeScript® RT reagent Kit (TaKaRa). Quantitative real-time PCR (qRT-PCR) was performed using ABI Stepone Plus instrument. Gene-specific primers and internal control primers were listed in Additional file 1: Table S6. All qRT-PCR reactions were performed in triplicates for each cDNA sample with an annealing temperature of 57°C and a total of 40 cycles of amplification. The relative expression levels were calculated by the comparative C T method.
Prediction of lncRNA function based on co-expression and genomic co-location A number of investigations have indicated that one major function of lncRNAs is regulating the expression of neighboring protein-coding genes via epigenetic modification or transcriptional co-activation/repression [13,43,44]. The relative loci between lncRNAs and their neighbors can be exhibited using Integrative Genomics Viewer [71]. Moreover, differentially expressed lncRNAs and mRNAs were forecasted to play roles in regulation of tolerance to osmotic and salt stresses. Therefore, the genomic colocational analysis of these lncRNAs and mRNAs was performed. We defined two genes as a co-expressed and co-located pair if they were co-expressed and spaced by less than 100 kb, according to the previously described method [72].
The neighbors of lncRNA genes were analyzed by Gene Ontology (GO) [73], and GO terms were enriched when significance (P) was less than 0.05 using Blast2GO [74]. Interaction networks among lncRNAs and proteincoding RNAs were constructed based on co-expression and genomic co-location using software Cytoscape [75].