Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Identification of tissue-specific and cold-responsive lncRNAs in Medicago truncatula by high-throughput RNA sequencing

BMC Plant Biology volume 20, Article number: 99 (2020) Cite this article

Abstract

Background

Long non-coding RNAs (lncRNAs) play important roles in the regulation of plant responses to environmental stress by acting as essential regulators of gene expression. However, whether and how lncRNAs are involved in cold acclimation-dependent freezing tolerance in plants remains largely unknown. Medicago truncatula is a prominent model for studies of legume genomics, and distinguished by its cold-acclimation characteristics. To determine the roles of lncRNAs in plant cold stress response, we conducted genome-wide high-throughput sequencing in the legume model plant M. truncatula.

Results

RNA-seq data were generated from twelve samples for the four treatments, i.e., non-cold treated leaves and roots, cold-treated leaves and roots of M. truncatula Jemalong A17 seedlings. A total of 1204 million raw reads were generated. Of them, 1150 million filtered reads after quality control (QC) were subjected to downstream analysis. A large number of 24,368 unique lncRNAs were identified from the twelve samples. Among these lncRNAs, 983 and 1288 were responsive to cold treatment in the leaves and roots, respectively. We further found that the intronic-lncRNAs were most sensitive to the cold treatment. The cold-responsive lncRNAs were unevenly distributed across the eight chromosomes in M. truncatula seedlings with obvious preferences for locations. Further analyses revealed that the cold-responsive lncRNAs differed between leaves and roots. The putative target genes of the lncRNAs were predicted to mainly involve the processes of protein translation, transport, metabolism and nucleic acid transcription. Furthermore, the networks of a tandem array of CBF/DREB1 genes that were reported to be located in a major freezing tolerance QTL region on chromosome 6 and their related lncRNAs were dissected based on their gene expression and chromosome location.

Conclusions

We identified a comprehensive set of lncRNAs that were responsive to cold treatment in M. truncatula seedlings, and discovered tissue-specific cold-responsive lncRNAs in leaves and roots. We further dissected potential regulatory networks of CBF Intergenic RNA (MtCIR1) and MtCBFs that play critical roles in response and adaptation of M. truncatula to cold stress.

Background

The discovery of non-coding RNAs (ncRNAs) including short (22–33 nucleotides) and long (> 200 nucleotides) ncRNAs has changed the traditional definition of a gene [1, 2]. Long non-coding RNAs (lncRNAs), which are distinguished by the lack of any obvious open reading frames (ORFs), are mainly transcribed by RNA Pol II, spliced, 5′-capped and even polyadenylated at 3′ end [2, 3]. In addition to RNA Pol II-derived lncRNAs, other classes of lncRNAs that were transcribed by two plant-specific DNA-dependent RNA polymerases, RNA Pol IV and RNA Pol V have also been reported [4]. LncRNAs can be classified into the following categories according to their genomic origins. (1) Antisense, when one or more exons of another transcript are overlapped on the opposite strand, respectively. (2) Intronic, when they are derived wholly from within an intron of a second transcript. (3) Intergenic, when they lie within the genomic intervals between two genes [2]. LncRNAs participate in the regulation of numerous biological phenomena, including those of imprinting genomic loci, shaping chromosome conformation and allosterically regulating enzymatic activity [3, 5]. LncRNAs function as key regulators of diverse mechanisms in biological processes, e.g., acting as scaffolds, decoys or signals through genomic targeting by cis or trans, and resulting in the down-regulation or overexpression of target genes [6, 7]. Studies on lncRNAs in plants have shown that they play important roles in a wide range of biological processes, especially in reproductive development and responses to environmental stresses [3, 8].

In plants, genome-wide identification of lncRNAs has been conducted in maize [9], Arabidopsis [10], Populus [11], Medicago truncatula [12], tomato [13] and others. The early functionally characterized plant lncRNAs are from the regulatory pathway of FLC, a master repressor gene involved in flowering in Arabidopsis. The lncRNAs of COLDAIR and COOLAIR negatively modulate FLC by different models [14,15,16]. A large number of lncRNAs involved in the regulation of plant responses to abiotic stresses has been characterized in recent years [17,18,19].

Plants grown in temperate and cold regions can enhance their tolerance to freezing by exposure to low, non-freezing temperatures for a certain period, referred to as cold acclimation [20, 21]. Numerous molecular changes during cold acclimation are responsible for cold acclimation-induced enhancement of freezing tolerance [22]. The CBF/DREB1 (C-repeat binding factor/dehydration-responsive element binding factor 1) activates the downstream CRT/DRE-containing cold-regulated (COR) genes, which in turn contributes to enhancing tolerance of plants to freezing [23,24,25]. ICE1 [26], ICE2 [27] and three closely related CAMTAs [28] have been identified as positive regulators of CBFs. The positive regulation of RDM4 on expression of AtCBF2 and AtCBF3 suggests that RDM4 is important for Pol II transcription of CBFs, and that it plays a critical role in tolerance of Arabidopsis to cold stress [29]. Recent studies have demonstrated that the cold-activated plasma membrane protein kinase CRPK1 phosphorylates 14–3-3 proteins, triggering its nuclear translocation to impair the stabilization of the transcription factor CBFs for a feedback of excessive cold defense response during cold stress in Arabidopsis [30]. The lncRNA SVALKA was identified as a negative regulator of CBF expression and plant freezing tolerance in Arabidopsis [31], however, how the CBF genes are activated by lncRNAs during cold acclimation remains to be explored.

Medicago truncatula is an annual forage crop [32] and has become a prominent model for legume genomics [33,34,35]. Given that M. truncatula is closely phylogenetically related to the common legume forage alfalfa (Medicago sativa), it is a valuable material to study molecular physiology of environmental stress in legume plants [36,37,38]. Alfalfa is a freezing tolerant legume species with great ability to cold acclimate, and capable of accumulating Cold-Acclimation-Specific (CAS) proteins during cold acclimation [36, 39, 40]. The CASs are homologous to COR genes that are cis-regulated by CBF/DREB1 factors [36, 38]. Similar to other plant species grown in temperate zones, M. truncatula plants have cold-acclimation traits [36, 41]. A tandem array of CBF/DREB1 genes was located in a major freezing tolerance QTL region on chromosome 6 of M. truncatula [42]. To identify cold-responsive lncRNAs in M. truncatula, we conducted genome-wide high-throughput sequencing for four treatments (Non-cold-treated leaves, Non-cold-treated roots, Cold-treated leaves, Cold-treated roots). We identified a comprehensive set of lncRNAs that were responsive to cold treatment in M. truncatula seedlings. Furthermore, a possible regulatory network of CBF Intergenic RNA (MtCIR1) and MtCBFs in M. truncatula was uncovered.

Results

Identification of cold-responsive lncRNAs in M. truncatula seedlings by high-throughput sequencing

To identify cold-responsive lncRNAs, we conducted RNA sequencing from 12 cDNA libraries with three repeats for four treatments (Non-cold-treated leaves, NT-leaves; Non-cold-treated roots, NT-roots; Cold-treated leaves, CT-leaves; Cold-treated roots, CT-roots) of M. truncatula seedlings. The cDNA libraries were constructed by synthetic adaptors using mRNAs isolated from leaves and roots of three-week-old M. truncatula seedlings by cold treatment at 4 °C and non-cold treatment at 26 °C for 5 h, respectively. Because of the low expression levels of lncRNAs in animals and plants [2], high depth RNA sequencing was performed to generate more than 1,204,140,634 raw reads from the 12 cDNA libraries (Table 1). To assess the quality of data acquired by RNA sequencing (RNA-seq), each base in the reads was assigned a quality score (Q) by a phred-like algorithm using FastQC [43]. The results showed that the data were highly credible with Q20 higher than 95% (Supplementary Table S1). After mapping sequencing results to the M. truncatula A17 genome (Supplementary Table S4), a comprehensive pipeline was constructed to identify unique lncRNAs. This led to identification of more than 30,000 mRNAs and 10,000 unique lncRNAs for each one of 12 cDNA libraries (Table 1; Additional file 8: Data Set S1; Additional file 9: Data Set S2; Additional file 10: Data Set S3; Additional file 11: Data Set S4).

Table 1 Statistical data of the RNA-Seq results for experimental samples
Full size table

To identify total numbers of lncRNAs in the four treatment samples, we merged data of three biological repeats (Correlation coefficient was shown in Supplementary Table S2 and S3) for each treatment regardless of repeatability, and 19,014, 16,298, 13,922 and 17,026 unique lncRNAs were identified in NT-leaves, NT-roots, CT-leaves and CT-roots, respectively (Fig. 1a). The numbers of identified lncRNAs in leaves were greater than those in roots of NT (19,014/16,298) scenarios, but they were less after cold treatment (13,922 vs 17,026) (Fig. 1a). The numbers of identified lncRNAs in CT-leaves were reduced by 5092, which resulted from 6494 disappearance and 1402 appearance compared with NT-leaves. In contrast, the numbers of identified lncRNAs were increased by 728 in CT-roots, which resulted from 2155 disappearance and 2883 appearance compared with NT-roots (Fig. 1a).

Fig. 1
figure 1

Identification of cold-responsive lncRNAs in M. truncatula seedlings. a: Numbers of common/specific lncRNAs identified in NT-L, CT-L, NT-R and CT-R. b: Total numbers of lncRNAs identified in leaves (NT, CT) and roots (NT, CT). Values shown in (a) and (b) were obtained from pooling data of three independent RNAseq experiments. Length distribution of lncRNAs in leaves (c) and roots (d). e: Numbers of lncRNAs containing different numbers of exon in NT-L, CT-L, NT-R and CT-R. f: Numbers of three type lncRNAs identified in NT-L, CT-L, NT-R and CT-R. Capital letters indicate a significant different at p < 0.05 according to t-test between three type lncRNAs in same treatment samples and small letters indicate a significant different at p < 0.05 according to t-test between NT-L, CT-L, NT-R and CT-R of same type lncRNAs. Values shown in (c), (d), (e) and (f) are means ± SE with three independent RNAseq experiments. Antisense-lncRNA: lncRNA overlapping with reference on the opposite strand; intronic-lncRNA: lncRNA falling entirely within a reference intron; incRNA: intergenic lncRNA; NT: non-cold treated; CT: cold treated; L: leaves; R: roots

Full size image

We identified 24,368 unique lncRNAs by pooling data from four treatments with 12 sample libraries (Fig. 1b). Among these lncRNAs, we simultaneously identified 10,934 (44.9%) unique lncRNAs in the four treatments, and 4465 (18.3%) and 3230 (13.3%) lncRNAs were specific to leaves and roots, respectively (Fig. 1b). In addition, 3208 (13.2%) and 516 (2.1%) lncRNAs were specifically identified in non-cold-treated and cold-treated leaves, and 824 (3.4%) and 1278 (5.2%) in non-cold treated and cold treated roots, respectively (Fig. 1b).

In terms of the length, the majority of lncRNAs was relatively short, such that the percentage of lncRNAs shorter than 2000 nt accounted for 81.1, 77.4, 80.5 and 80.8% in NT-leaves, NT-roots, CT-leaves and CT-roots, respectively (Fig. 1c and d). Specifically, we found that lncRNAs shorter than 400 nt were dominant ones. For example, lncRNAs shorter than 400 nt were 42.6, 31.6, 40.8 and 41.8% of total lncRNAs in NT-leaves, NT-roots, CT-leaves and CT-roots, respectively (Fig. 1c and d).

We also analyzed exon numbers in the lncRNAs. As shown in Fig. 1e, the numbers of lncRNAs dramatically reduced with increasing exon numbers, and the numbers of lncRNAs containing only one exon were the most, and accounted for about 42.4, 39.8, 31.5 and 41.0% in NT-leaves, NT-roots, CT-leaves and CT-roots, respectively (Fig. 1e).

In the present study, we classified the lncRNAs into three categories, i.e., antisense-lncRNA (lncRNA overlapping with reference on the opposite strand), intronic-lncRNA (lncRNA falling entirely within a reference intron) and lincRNA (intergenic lncRNA) according to their genomic origins (Fig. 1f). Among the three types of lncRNAs, the most and least numbers belonged to intronic-lncRNA and antisense-lncRNA ones in both leaves and roots (Fig. 1f).

Analyses of lncRNAs distribution on chromosome revealed that the identified lncRNAs were unevenly distributed across the eight chromosomes in M. truncatula seedlings with obvious preferences for locations (Additional file 2: Fig. S1; Additional file 3: Fig. S2).

Characterization of cold-responsive lncRNAs in M. truncatula seedlings

We characterized the cold-responsive lncRNAs according to the following five aspects: change fold (FC(log2)), exon number, length of lncRNAs, distribution on chromosome, and type of lncRNAs in leaves and roots of M. truncatula. The lncRNAs displaying different expression were selected based on change fold≥2 of TPM (Transcripts Per Million) and p-value< 0.05 for three biological repeats between cold and non-cold treated leaves and roots (Additional file 8: Data Set S1).

Among identified lncRNAs, 983 and 1288 were responsive to cold treatment in the leaves and roots, respectively (Fig. 2a). The numbers of cold-responsive lncRNAs were greater in roots than those in leaves as well as those of up- or down-regulated lncRNAs by cold treatment (Fig. 2a). The numbers of lncRNAs that were up-regulated were greater than those of down-regulated lncRNAs in both leaves and roots (Fig. 2b). Moreover, only a few lncRNAs were specifically found to be up-regulated by cold treatment in cold treated leaves and roots (Additional file 8: Data Set S1). As shown in Fig. 2b, the numbers of lncRNAs with expression levels of up- or down-regulated were reduced with increasing (FC(log2)) in both leaves and roots. The highest numbers of up- or down-regulated lncRNAs were those whose expression changes were FC(log2) < 4. For example, they were 91.6 and 89.2% for up-regulated and down-regulated in leaves, while they were 93.6 and 95.2% for up-regulated and down-regulated in roots, respectively (Fig. 2b).

Fig. 2
figure 2

Characteristics of cold-responsive lncRNAs in M. truncatula seedlings. a: Numbers of cold-responsive lncRNAs in leaves and roots of M. truncatula seedlings. b: Numbers of lncRNAs with different FC(log2) which were regulated by cold treatment in leaves and roots of M. truncatula seedlings. Data on X axis indicated varying (FC(log2)) values ranging from ≥1 to ≥10. c: Numbers of lncRNAs containing different exon numbers in response to cold treatment in leaves and roots of M. truncatula seedlings. d: Numbers of lncRNAs with different length responded to cold treatment in leaves and roots of M. truncatula seedlings. e: Numbers of lncRNAs distributed on chromosomes in response to cold treatment in leaves and roots of M. truncatula seedlings. f: Densities of lncRNAs on 8 chromosomes responded to cold treatment in leaves and roots of M. truncatula seedlings. CT-res: cold treatment responsive-lncRNAs; up- or down-regulated: expression of lncRNAs induced or reduced by cold treatment; L: leaves; R: roots. Three types of lncRNAs (antisense-lncRNA, intronic-lncRNA and lincRNA) regulated by cold treatment in leaves (g) and roots (h) of M. truncatula seedlings. Antisense-lncRNA: lncRNA overlapping with reference on the opposite strand; intronic-lncRNA: lncRNA falling entirely within a reference intron; lincRNA: intergenic lncRNA

Full size image

The numbers of lncRNAs showing up- or down-regulation were also reduced with increasing exon numbers except lncRNAs containing two exons in either leaves or roots (Fig. 2c). Different from the results that the numbers of identified lncRNAs containing only one exon were most in the four treatments, the most numbers of cold-responsive lncRNAs were those containing two exons (Fig. 2c). Although the majority (about 40%) of lncRNAs identified in M. truncatula seedlings was shorter than 400 nt (Fig. 2c and d), lncRNAs shorter than 2000 nt showed comparable responses to cold treatment regardless of their length (Fig. 2d).

The numbers of up-regulated lncRNAs were more than those of down-regulated lncRNAs in leaves across the eight chromosomes in response to cold treatment (Fig. 2e and f). In contrast, the cold-responsive lncRNAs exhibited different patterns in roots. For example, the numbers of up-regulated lncRNAs in roots were less than those of down-regulated lncRNAs on chr 6 and chr 8, while the numbers of up-regulated lncRNAs were more than those of down-regulated lncRNAs on the remaining chromosomes in response to cold treatment (Fig. 2e and f). Further analyses indicated that cold treatment-induced increases in number and density of identified lncRNAs differed among 8 chromosomes in leaves and roots (Fig. 2e and f). The largest number was found on chr3 (17.0%), while the fewest one was on chr6 (6.4%) in leaves (Fig. 2e). The most number of lncRNAs up-regulated by cold treatment was on chr4 (20.2%), while the fewest one was on chr6 (21.1%) in roots (Fig. 2e). Although the chromosomes with the highest density of up- and down-regulated lncRNAs in leaves differed from those in roots, the chromosomes with the lowest density of up- and down- regulated lncRNAs were on chr3 in both leaves and roots (Fig. 2e and f). These results suggest that each chromosome has distinct patterns in response to cold treatment in term of expression of lncRANs.

The most numbers of cold-responsive lncRNAs belonged to intronic-lncRNAs both in leaves and roots, but the least ones were linc-lncRNAs in leaves and antisense-lncRNAs in roots (Fig. 2g and h). Moreover, different types of lncRNAs differed in their responses to cold treatment, such that the expression of 49.4 and 45.0% of intronic-lncRNAs in leaves and roots was up-regulated by cold treatment, but it was only 6.5 and 4.9% for antisense-lncRNAs, and 4.5 and 4.9% for linc-lncRNAs, respectively (Fig. 2g and h). Similar to the up-regulated lncRNAs, the numbers of down-regulated intronic-lncRNAs were also most, but the least ones were antisense-lncRNAs in both of leaves and roots (Fig. 2g and h). These results may suggest that intronic-lncRNAs are most sensitive lncRNAs to cold treatment.

Functional characterization of cold-responsive lncRNAs in M. truncatula seedlings

To uncover potential functions of the cold-responsive lncRNAs, we analyzed Gene Ontology (GO) terms of genes that were co-expressed and co-localized with the cold-responsive lncRNAs based on the transcriptional expression of three biological repeats between cold and non-cold treated leaves and roots and genomic location in 100 kb. The top 10 terms of biological processes, cellular components and molecular functions were analyzed according to significant enrichments (corrected P < 0.05) in leaves and roots (Fig. 3; Additional file 1: Table S5 and S6).

Fig. 3
figure 3

Gene Ontology (GO) analysis. GO terms of genes that were co-expressed and co-localized with the cold-responsive lncRNAs based on transcriptional expression and genomic location in leaves (a) and roots (b) of M. truncatula seedlings. The reliability was calculated by –log10 (p < 0.05). Black color words in Y-axis indicated the genes of these GO terms were up- or down-regulated by cold treatment. Orange color words in Y-axis indicated the genes of these GO terms were up-regulated by cold treatment. Blue color words in Y-axis indicated the genes of these GO terms were down-regulated by cold treatment. More information is detailed in Additional file 1: Table S5 and S6

Full size image

Among the top 10 terms of biological processes, expression of all genes associated with DNA-templated transcription (initiation) (GO:0006352) in leaves and those associated with base-excision repair (GO:0006284) in roots was induced by cold treatment. Expression of all genes involved in base-excision repair (GO:0006284) and vesicle-mediated transport (GO:0007018) in leaves and genes in nucleotide-excision repair (GO:0006289) in roots was reduced by cold treatment. Genes in the majority of GO groups were up- or down- regulated by cold treatment (Fig. 3). For the top 10 terms of cellular components in leaves, cold treatment up-regulated expression of all genes in the three GO groups of leaves, i.e., ribosome, small-subunit processome and signal peptidase complex (GO:0005840, GO:0032030, GO:0005787). In contrast, cold treatment led to down-regulation of expression of all genes in nucleosome (GO:0000786) and mitochondrial outer membrane (GO:0005741) in roots. Genes in other GO groups of cellular component were up- or down regulated by cold treatment in either leaves or roots (Fig. 3). The molecular function of genes co-expressed and co-localized with cold-responsive lncRNAs in leaves and roots were extremely complicated. Expression of genes in galactosyl transferase activity group (GO:0008378) was induced, and expression of genes in microtubule motor activity group (GO:0003777) was reduced by cold treatment (Fig. 3). Proteins encoded by these mRNAs had different functions in leaves and roots, suggesting organ-specific responses to cold treatment. These results highlight the diversity of lncRNAs genomic location and complexity of regulation function acting on their potential target genes.

To further characterize the expression relationship between lncRNAs and their potential target protein-coding genes, correlations between expression patterns of lncRNAs and their potential target genes that were responsive to cold treatment were analyzed (Fig. 4). As shown in Fig. 4, more than 30% of potential target genes displayed similar expression patterns with intronic-lncRNAs and antisense-lncRNAs (Fig. 4a and b). For example, about 51.7% potential target genes of cold-induced up-expression of intronic-lncRNAs in leaves (L-up) were detected, while 43.6% potential target genes of cold-induced expression of intronic-lncRNAs in roots (R-up) were found. Similar to cold-induced expression of lncRNAs, 31.3% potential target genes of L-down (i.e., expression was reduced by cold treatment in leaves) intronic-lncRNAs were reduced by cold treatment, while it was 43.1% for potential target genes of R-down (i.e., expression was reduced by cold treatment in roots) intronic-lncRNAs (Fig. 4a). Similar scenarios were found for those potential target genes of antisense-lncRNAs (Fig. 4b). In addition, the expression patterns of potential target genes related with lincRNAs were more complex (Fig. 4c and d). The expression patterns of upstream potential target genes exhibited little similarity to their related lincRNAs (Fig. 4c). On the contrary, the expression of downstream potential target genes was positively related with expression of adjacent lincRNAs (Fig. 4d). For example, about 40.8% downstream potential target genes of L-up lincRNAs were induced by cold treatment, while 38.8% downstream potential target genes of R-up lincRNAs were induced by the same cold treatment. Similar to cold-induced expression of lncRNAs, 43.8% downstream potential target genes of down-regulated lincRNAs in leaves (L-down) were reduced by cold treatment, while 35.1% downstream potential target genes of down-regulated lincRNAs in roots (R-down) (Fig. 4d).

Fig. 4
figure 4

Correlation analyses of lncRNAs and their potential target genes in M. truncatula. Expression models of intronic-lncRNAs (a), antisense-lncRNAs (b) and their potential target genes in leaves and roots of M. truncatula seedlings responsive to cold treatment. Expression models of lincRNAs and their potential upstream target genes (c), downstream genes (d) in leaves and roots of M. truncatula seedlings responding to cold treatment. The expressions were analyzed by RNA-seq (TPM). Antisense-lncRNA: lncRNA overlapping with reference on the opposite strand; intronic-lncRNA: lncRNA falling entirely within a reference intron; lincRNA: intergenic lncRNA; L: leaves; R: roots; up: expression were induced by cold treatment; down: expression were reduced by cold treatment; ND: no detected

Full size image

Analyses of lncRNA-MtCBFs networks

C-repeat/DRE binding proteins (CBFs) have been identified to be pivotal transcription regulatory factors of cold-responsive (COR) genes during cold treatment in many plant species [22]. A tandem array of CBF/DREB1 genes has been reported to be located in a major freezing tolerance QTL region on chromosome 6 in M. truncatula [42]. To determine the roles of cold-responsive lncRNAs in regulating MtCBFs, we monitored the expression patterns of lncRNAs and MtCBFs by taking into account of their locations on the chromosomes. One reverse-direction transcription lncRNA (LNC_016398-MtCIR1, M.truncatulaCBFs Intergenic RNA) (Fig. 5a and c) and 7 CBF genes (Mt6g465420, Mt6g465430, Mt6g465450, Mt6g465460, Mt6g465510, Mt6g465530 and Mt6g465690) (Fig. 5b and d) were induced by cold treatment by RNA-seq (Fig. 5a and b), and verified by real-time quantitative Q-PCR (Fig. 5c and d). Analyses of genomic location of the MtCIR1 and seven CBF genes indicated that they were neighborly distributed on chromosome 6 (Fig. 5e). MtCIR1 (733 bp) was located in the intergenic region of Mt6g465420 and Mt6g465430, and it did not show any overlap with these genes (Fig. 5e).

Fig. 5
figure 5

Relationship between lncRNA MtCIR1 and their potential target MtCBF genes in M. truncatula plants. Transcriptional responses of MtCIR1 (a) and MtCBF genes (b) in M. truncatula plants to cold treatment at 4 °C for 5 h analyzed by RNA-seq (TPM). Transcription responses of MtCIR1 (c) and MtCBF genes (d) in M. truncatula plants to cold treatment at 4 °C for 5 h verified by real-time quantitative Q-PCR. NT: non-cold treated; CT: cold treated; L: leaves; R: roots. Values are means ± SE (n = 3) and the asterisks indicate significant differences at p < 0.05 between cold-treated sample (CT) and non-cold-treated control (NT) according to t-test. e: Genome map of the MtCBF gene clusters and MtCIR1. MtCIR1 and MtCBF gene clusters are located in a 145.3 kb on chromosome 6 of M. truncatula Jemalong A17. Red rectangle indicates exon. Red arrow indicates exon and direction of gene from 5′ to 3′. Green line indicates intron. Blue rectangle indicates lncRNA. Black line indicates intergenic region. Black numbers indicate length (bp) of intergenic region

Full size image

lncRNA MtCIR1 was not only closely neighbored by the seven MtCBFs on the chromosome 6 (Fig. 5e), its expression level was also comparable to that cold-induced MtCBFs (Fig. 5a-d). To confirm the results, we tested transcript responses of MtCIR1 and their potential target MtCBF genes to cold treatment by quantitative real-time PCR in leaves of M. truncatula (Fig. 6). The expression of MtCIR1 was up-regulated by cold treatment (Fig. 6a). For instance, after cold treatment for 5 h, transcripts of MtCIR1 was increased by 7 fold (Fig. 6a). Transcripts of the seven MtCBFs that were located on chr6 and neighbored by the MtCIR1 were transiently increased in leaves (Fig. 6b). An important question is whether the cold treatment-induced expression of MtCIR1 preceded the expression of MtCBFs. To answer this question, we monitored the expression of MtCIR1 and MtCBFs during early stage (0 to 2 h) of cold treatment in leaves. As shown in Fig. 6a, the transcript levels of MtCIR1 were induced markedly within 2 h of cold treatment. The cold-induced increase in MtCIR1 was followed at 5 h by accumulation of transcripts for MtCBFs (Fig. 6b). In contrast to other four MtCBFs, which peaked at 5 h and then fall during cold treatment, the expression of three MtCBFs, Mt6g465510, Mt6g465530 and Mt6g465690 also peaked at 5 h and remained at the high level during cold treatment up to 24 h (Fig. 6b). The results that up-expression of MtCIR1 was followed by the induction of MtCBF transcripts during cold treatment may suggest the existence of crosstalk between the lncRNA MtCIR1 and MtCBF genes.

Fig. 6
figure 6

Responses of lncRNA MtCIR1 and their potential target MtCBF genes to cold treatment. a: Transcriptional responses of MtCIR1 in leaves of M. truncatula seedlings subjected to cold treatment. b: Transcriptional responses of MtCBF genes in leaves of M. truncatula seedlings subjected to cold treatment. Leaf RNA was extracted from 3-week-old seedlings cold treated at 4 °C for varying time periods and then the relative expression was detected using real-time quantitative RT-PCR. Values are means ± SE (n = 3) and the asterisks indicate significant differences at p < 0.05 according to t-test relative to room temperature control (0 h)

Full size image

Discussion

Studies on lncRNAs in plants have shown that they play important roles in a wide range of biological processes, especially in regulating plant responses to biotic and abiotic stress, such as drought stress in maize [9] and Populus trichocarpa [11], salt and drought stress in Medicago truncatula [12], P. infestans-resistant in tomatoes [13]. However, we know little about whether involvement of lncRNAs in cold acclimation-dependent freezing tolerance in legume model plant Medicago truncatula. In the present study, we identified 24,368 unique lncRNAs that is similar to numbers obtained in M. truncatula induced by salt and drought stress lncRNAs [12]. Among 24,368 lncRNAs, 983 and 1288 were responsive to cold treatment in the leaves and roots. The cold-responsive lncRNAs found in our study are more than those identified in cassava, and this may result from the different sequencing methods and plant species, such as different genome sizes and sensitivities to low temperature stress [44].

Our results showed that total numbers of identified unique lncRNAs and mRNAs in leaves were greater than those in roots, but the proportion of cold-responsive lncRNAs and mRNAs was higher in roots than that in leaves (Figs. 1 and 2; Additional file 4: Fig. S3; Additional file 5: Fig. S4). These results may suggest that different organs have specific responses to cold treatment, such that lncRNAs and mRNAs in roots may be more sensitive to cold treatment. The higher proportion of cold-responsive lncRNAs and mRNAs in roots relative to those in leaves may imply the greater importance of roots in response to cold treatment. However, few studies have focused on physiological and molecular responses of roots to cold stress so far [45, 46]. It has been suggested that the capacity for cold acclimation and frost resistance in four annual legumes were related to a higher root to shoot ratio and higher concentrations of solutes in roots [47]. Therefore, root-mediated physiological processes may play important roles in the regulation of cold treatment-dependent freezing tolerance in plants.

Our results that the majority of GO terms was involved in the regulation of varying biological processes are consistent with the reports suggesting existence of complicated mechanisms in plant cold stress response [21, 48]. GO terms involved in defense response may directly participate in the protection of cells under stress conditions. The GO term of 0006355 associated with the regulation of transcription has been reported to play a role in response to freezing stress in M. sativa [49]. It has been reported that the gene coding for Cu/Zn superoxide dismutase, which is cleaved by miR398 in plants, was involved in the regulation of cold-stress responses by acting as a scavenger of reactive oxygen species (ROS) [50, 51]. These results may imply that correlations between lncRNAs and microRNAs may exist in the modulation of responses to cold stress. An interesting finding is that genes included in base-excision repair (GO: 0006284) were down-regulated in leaves, but up-regulated in roots by cold treatment (Fig. 3). DNA damaging agents can be a great threat to genomes. To protect DNA from the damage, a number of DNA repairing strategies have been developed. The base- excision-repair pathway is a highly conserved mechanism in the organisms, and has been suggested to account for removal and repair of mutagenic oxidative DNA lesions [52]. Although there are extensive studies on the roles of base excision repair system in stress and plant hormone signaling [53, 54], little is known about their effects on plant cold stress response. The GO terms of base-excision repair identified in this research indicate that the base-excision repair may play important roles in repairing DNA during cold treatment.

As key transcriptional activators, the gene family of dehydration and cold response (CBF/DREB1) activates the downstream cold-regulated (COR) genes, which in turn contributes to enhancement of tolerance to freezing stress [22, 24, 55]. A major freezing tolerance QTL (Mt-FTQTL6) accounting for 40% of the phenotypic variation among 15 M. truncatula accessions has been mapped to a region of M. truncatula chromosome 6 [56]. A tandem array of CBF/DREB1 genes was located in the major freezing tolerance QTL (Mt-FTQTL6) region on the chromosome 6 in M. truncatula [42]. However, as an early event of a low temperature-stimulated signaling cascade, how the expression of CBFs is induced by cold treatment remains largely unclear. Numerous reports have demonstrated that lncRNAs participate in the regulation of various biological processes by interacting with DNA and RNA molecules, and transcription factors, leading to alterations of target genes [5]. Furthermore, the expression patterns of lncRNAs are often correlated with those of mRNA in both cis and trans manners, suggesting that certain lncRNAs may be co-regulated in the expression networks [57]. In the present study, we found the involvement of lncRNAs and MtCBFs genes (GO: 0005634) in cold response by RNA-sequencing and gene expression analyses. We demonstrated that expression of six MtCBFs was related to the expression of one lncRNA MtCIR1. Specifically, we found that MtCIR1 was an intergenic noncoding RNA located in a close proximity to MtCBF genes, and that the expression of the MtCIR1 and MtCBF genes was induced by cold treatment. This observation is in line with a report demonstrating that about half of the intergenic noncoding RNAs is transcribed close to protein-coding genes [57]. Therefore, these data may provide important clues for further dissection about the transcription regulation of intergenic noncoding RNAs on neighboring genes during cold treatment.

The observations that transcript levels of lncRNA MtCIR1 increased within 2 h of exposure to low temperature, followed by accumulation of MtCBFs at 5 h may suggest a regulatory network between MtCBFs and MtCIR1. However, these results differ from previous studies that the accumulation of AtCBF1–3 transcripts was detected within 15 min of plants upon exposure to low temperature, followed by accumulation of COR gene transcripts at 2 h [58]. Despite the high similarities in sequence and close evolution relationships between the six MtCBFs identified in this research with AtCBF1–3, the difference in cold-response may imply that the sensitivity of CBFs transcript to cold treatment is specific to plant species. For example, Arabidopsis plants are more sensitive to cold stress than M. truncatula plants (Additional file 6: Fig. S5; Additional file 7: Fig. S6). Recent studies found an lncRNA SVALKA in a cold-sensitive region of the Arabidopsis genome. Mutations in SVALKA affect CBF1 expression and freezing tolerance [31]. Whether lncRNA MtCIR1 in M. truncatula has similar function in the regulation of plant freezing tolerance with SVALKA in Arabidopsis warrants further studies using genetics and physiology methods.

Many studies have found that ICE1, ICE2, the three closely related CAMTAs, and RDM4, are positive regulators of CBFs, and that the phosphorylated 14–3-3 proteins destabilize CBF proteins [26,27,28,29,30, 59]. However, whether MtCIR1 is correlated with these proteins in the regulation of MtCBFs and freezing tolerance remains to be elucidated. Shu et al. [49] identified a number of freezing- and Medicago-specific miRNAs involved in the regulation of freezing tolerance in M. sativa. One important mechanism of lncRNAs at the post-transcriptional is to function as mediators, including functioning as the precursors of small RNAs, and acting with miRNAs to regulate mRNA turnover [60]. Whether the cold responsive MtCBFs and lncRNA MtCIR1 found in the present study is also related with these microRNA remains to be uncovered in the future study.

The expression levels of sequencing and qPCR were inconsistent in our research, especially at expression levels of lncRNAs (Fig. 5a-d). There are mainly three factors that may account for these results. Firstly, a certain degree of inconsistency between the results obtained by the two methods may be expected, especially for the low expression of RNAs because the sequencing and qPCR are two different detection methods. Secondly, qPCR was used to verify the sequencing results, which verified the difference trend, i.e. up-regulated or down-regulated, rather than the difference multiple. Lastly, duplication during sequencing might partly interfere quantitative results of RNA transcripts [61, 62].

Conclusions

In the present study, we identified a large number of cold-responsive lncRNAs in both leaves and roots of legume model plant M. truncatula seedlings. We further demonstrated that the cold-responsive lncRNAs were tissue-specific, i.e., the numbers of identified lncRNAs in leaves were greater than those in roots of non-cold stress scenarios, but they were less after cold treatment. We found that the three types of intronic-lncRNAs differed in their responses to cold stress, with the intronic-lncRNAs being the most sensitive to cold stress. Another interesting finding is that cold-responsive lncRNAs were unevenly distributed across the eight chromosomes in the genome of M. truncatula with obvious preferences for locations. Furthermore, we dissected a regulatory network of lncRNA-MtCBFs that may play a critical role in response and adaptation of plants to cold stress by integrating lncRNA MtCIR1 with their potential target MtCBF genes.

Methods

Plant materials and growth conditions

Seeds of M. truncatula, Jemalong A17 used in this work were kindly provided by Dr. Carroll Vance, USDA-ARS, Plant Sceience Research, St. Paul, MN, USA. The Jemalong A17 is the model plant of legume plants, whose genome has been sequenced [35, 63]. Seeds of Medicago truncatula were treated with concentrated sulfuric acid for 4 min, and then thoroughly rinsed with water. After chilled at 4 °C for 2 d, seeds were sown on 0.8% agar to germinate at 25 °C till the radicals being about 2 cm. The seedlings were planted in the plastic buckets (three seedlings in one bucket) filled with aerated nutrient solution. The composition of full-strength nutrient solution is: 2.5 mM KNO3, 0.5 mM KH2PO4, 0.25 mM CaCl2, 1 mM MgSO4, 100 μM Fe-Na-EDTA, 30 μM H3BO3, 5 μM MnSO4, 1 μM ZnSO4, 1 μM CuSO4 and 0.7 μM Na2MoO4 with pH of 6.0. Plants were grown in green house under 26 °C day/22 °C night, and 14-h photoperiod, 120 μmolm− 2 s− 1 conditions.

Cold treatment and sample collection

Nine three-week-old Jemalong A17 seedlings grown in three plastic buckets were cold treated at 4 °C, 60 μmolm− 2 s− 1 for 5 h (CT). In the meanwhile, nine three-week-old Jemalong A17 seedlings grown in another three plastic buckets transferred to 26 °C, 60 μmolm− 2 s− 1 for 5 h were used as control (NT). Nine leaflets (randomization approach) for leaf-samples and three entire roots for root-samples from three plants grown in one plastic bucket were collected for one biological replicate, respectively. Three biological replicates were maintained for sample collection by a temperature treatment was repeated three separate times in a single growth chamber under same condition. Finally, twelve samples for four treatments of non-cold treated leaves and roots (NT-L, NT-R), cold-treated leaves and roots (CT-L, CT-R) were collected for RNA sequencing. In the meanwhile, twelve sub-samples were collected for qRT-PCR test. After harvest, samples were frozen in liquid nitrogen, and stored at − 80 °C for RNA extraction.

RNA extraction and cDNA library construction

Total RNAs were extracted from different samples using the Trizol (Invitrogen) according to the manufacturer’s protocols. After DNA digestion with RNase-free DNase I (Promega). RNA degradation and contamination were evaluated on 1% agarose gels. The purity of RNA was checked by the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was determined by Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was checked by the RNA Nano 6000 Assay Kit (Bioanalyzer 2100, Agilent Technologies, CA, USA).

For each sample, an amount of 3 μg RNA was used for the preparation of RNA samples. Ribosomal RNA was first removed by Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA), and rRNA free residue was cleaned by ethanol precipitation. Thereafter, the sequencing libraries were generated using the rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the protocols described by the manufacturer. Products were purified (AMPure XP system) and library quality was determined with the Agilent Bioanalyzer 2100 system.

Sequencing and data analysis

The libraries were sequenced on an Illumina Hiseq 4000 platform and 150 bp paired-end reads were generated. The raw reads from the 12 samples were used for quality filtering. Clean reads were obtained by removal of reads containing adapter, ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC contents of the clean data were determined [43]. The low quality reads (Phred score < 20; read length < 50 bases) and reads with adapter contamination were removed to generate a set of high quality reads termed as clean data thereafter. All the downstream analyses were based on the clean data with high quality.

Reads mapping and transcriptome assembling

The clean reads mapped to the M. truncatula genome. Index of the Medicago Genome Sequences V4.0 (http://www.medicagohapmap.org/tools/blastform) was built using bowtie2 v2.2.8 and paired-end clean reads were aligned to the M. truncatula genome using HISAT2 v2.0.4 [64]. To construct transcriptome, the mapped reads were assembled de novo using Cufflinks [65]. All transcripts were required to be > 200 bp in length.

The mapped reads of each sample were assembled by StringTie (v1.3.1) [66] in a reference-based approach. StringTie uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus.

Coding potential analysis and identification of lncRNAs

Predication of transcripts with coding potential was made using the following tools, including CNCI (Coding-Non-Coding-Index) (v2)) [67], CPC (Coding Potential Calculator) (0.9-r2) [68], Pfam Scan (v1.3) [69] and PhyloCSF (phylogenetic codon substitution frequency) (v20121028) [70]. The left unknown transcripts without coding potential were taken as candidate set of lncRNAs.

Target gene prediction

We first identified coding genes that were located 100 k upstream and downstream of lncRNA, and then analyzed their functions. The genes from different samples were clustered with WGCNA (Weighted Gene Co-expression Network Analysis) [71] to search for common expression modules, and their function was analyzed by the functional enrichment analysis.

Quantification of gene expression level and analysis of differential expression

Kallisto-sleuth pipelines was used to calculate TPMs of both lncRNAs and coding genes in each sample [72]. TPM means transcript per kilobase of exon per million fragments mapped, calculated based on the length of the transcript and reads count mapped to this transcript. The differentially expressed genes (DEGs) with log2 fold change ≥2 (induced) and/ or ≤ − 2 (reduced) and a P value less than 0.05 for either of the sample in each pair wise comparison were considered to be significantly differentially expressed.

Quantitative real-time PCR (qRT-PCR)

Total RNAs were extracted using the Trizol (Invitrogen) following the protocols provided by the manufacturer. The DNA was digested with RNase-free DNase I (Promega). RNA of about 0.5 μg was reverse-transcribed into first-strand cDNA with PrimeScript® RT reagent Kit (TaKaRa). Quantitative real-time PCR (qRT-PCR) was conducted by ABI Stepone Plus instrument. Gene-specific and internal control primers were given in Supplementary Table S7. We performed three independent experiments from three biological repeats for qRT-PCR, and three measurements were made for each cDNA with an annealing temperature of 56 °C and a total of 40 cycles of amplification. The relative expression levels were calculated by the comparative CT method.

Gene ontology (GO) enrichment analysis

Analyses of Gene Ontology (GO) enrichment for differentially expressed genes or lncRNA target genes were performed by the GOseq R package, with correction of gene length bias [73]. GO terms with corrected P values less than 0.05 were taken as the significantly enriched by differentially expressed genes. We constructed interaction-networks among lncRNAs and protein-coding RNAs based on co-expression and genomic co-location by the software Cytoscape [74].

Statistical analysis

All experiments in this study were repeated independently at least three times. The results are given means ± SE. The statistical analysis was performed using SPSS17.0 software (Chicago, IL, USA). The t-test was used to determine whether effects of treatments were statistically different at p < 0.05 level.

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files.

Abbreviations

Antisense-lncRNA:

LncRNA overlapping with reference on the opposite strand

CAMTAs:

Calmodulin binding Transcription Activators

CAS:

Cold-Acclimation-Specific

CBF/DREB1 :

C-Repeat Binding Factor/Dehydration-Responsive Element Binding factor 1

CNCI:

Coding-Non-Coding-Index

COR:

CRT/DRE-containing Cold-regulated

CPC:

Coding Potential Calculator

CT:

Cold Treated

DEG:

Differentially expressed gene

FC(log2):

Log2-transformed fold change

GO:

Gene Ontology

IncRNA:

Intergenic lncRNA

Intronic-lncRNA:

lncRNA falling entirely within a reference intron

lncRNAs:

Long noncoding RNAs

MtCIR1 :

M. truncatula CBF Intergenic RNA 1

NT:

Non-cold Treated

ORFs:

Open Reading Frames

PhyloCSF:

Phylogenetic Codon Substitution Frequency

QC:

Quality Control

RDM4 :

RNA-directed DNA methylation 4

RNA-seq:

RNA sequencing

TPM:

Transcript per kilobase of exon per million fragments mapped

WGCNA:

Weighted Gene Co-expression Network Analysis

References

  1. Mattick JS, Makunin IV. Non-coding RNA. Hum Mol Genet. 2006;15:R17–29.

    Article  CAS  PubMed  Google Scholar 

  2. Ponting CP, Oliver PL, Reik W. Evolution and functions of long noncoding RNAs. Cell. 2009;136:629–41.

    Article  CAS  PubMed  Google Scholar 

  3. Chekanova JA. Long non-coding RNAs and their functions in plants. Curr Opin Plant Biol. 2015;27:207–16.

    Article  CAS  PubMed  Google Scholar 

  4. Wierzbicki AT, Haag JR, Pikaard CS. Noncoding transcription by RNA polymerase pol IVb/pol V mediates transcriptional silencing of overlapping and adjacent genes. Cell. 2008;135:635–48.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Chen LL, Carmichael GG. Decoding the function of nuclear long non-coding RNAs. Curr Opin Cell Biol. 2010;22:357–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Kornienko AE, Guenzl PM, Barlow DP, Pauler FM. Gene regulation by the act of long non-coding RNA transcription. BMC Biol. 2013;11:59.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Quinn JJ, Chang HY. Unique features of long non-coding RNA biogenesis and function. Nat Rev Genet. 2016;17:47–62.

    Article  CAS  PubMed  Google Scholar 

  8. Nejat N, Mantri N. Emerging roles of long non-coding RNAs in plant response to biotic and abiotic stresses. Crit Rev Bio technol. 2018;38:93–105.

    Article  CAS  Google Scholar 

  9. Zhang W, Han ZX, Guo QL, Liu Y, Zheng YX, Wu FL, et al. Identification of maize long non-coding RNAs responsive to drought stress. PLoS One. 2014;9:e98958.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  10. Liu J, Jung C, Xu J, Wang H, Deng SL, Bernad L, et al. Genome-wide analysis uncovers regulation of long intergenic noncoding RNAs in Arabidopsis. Plant Cell. 2012;24:4333–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Shuai P, Liang D, Tang S, Zhang ZJ, Ye CY, Su YY, et al. Genome-wide identification and functional prediction of novel and drought-responsive lincRNAs in Populus trichocarpa. J Exp Bot. 2014;65:4975–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Wang TZ, Liu M, Zhao MG, Chen RJ, Zhang WH. Identification and characterization of long non-coding RNAs involved in osmotic and salt stress in Medicago truncatula using genome-wide high-throughput sequencing. BMC Plant Biol. 2015;15:131.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Cui J, Luan YS, Jiang N, Bao H, Meng J. Comparative transcriptome analysis between resistant and susceptible tomato allows the identification of lncRNA16397 conferring resistance to Phytophthora infestans by co-expressing glutaredoxin. Plant J. 2017;89:577–89.

    Article  CAS  PubMed  Google Scholar 

  14. Liu F, Marquardt S, Lister C, Swiezewski S, Dean C. Targeted 30 processing of antisense transcripts triggers Arabidopsis FLC chromatin silencing. Science. 2009;327:94–7.

    Article  PubMed  CAS  Google Scholar 

  15. Heo JB, Sung S. Vernalization-mediated epigenetic silencing by a long intronic noncoding RNA. Science. 2011;331:76–9.

    Article  CAS  PubMed  Google Scholar 

  16. Castaings L, Bergonzi S, Albani MC, Kemi U, Savolainen O, Coupland G. Evolutionary conservation of cold-induced antisense RNAs of FLOWERING LOCUS C in Arabidopsis thaliana perennial relatives. Nat Commun. 2014;5:4457.

    Article  CAS  PubMed  Google Scholar 

  17. Qin T, Zhao HY, Cui P, Albesher N, Xiong LM. A nucleus-nocalized long non-coding RNA enhances drought and salt stress tolerance. Plant Physiol. 2017;175:1321–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Wang TZ, Zhao MG, Zhang XX, Liu M, Yang CG, Chen YH, et al. Novel phosphate deficiency-responsive long non-coding RNAs in the legume model plant Medicago truncatula. J Exp Bot. 2017;68:5937–48.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Wang Y, Luo XJ, Sun F, Hu JH, Zha XJ, Su W, et al. Overexpressing lncRNA LAIR increases grain yield and regulates neighbouring gene cluster expression in rice. Nat Commun. 2018;9:3516.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  20. Graham D, Patterson BD. Responses of plants to low, non freezing temperatures: proteins, metabolism, and acclimation. Annu Rev Plant Physiol. 1982;33:347–72.

    Article  CAS  Google Scholar 

  21. Guy CL. Cold acclimation and freezing stress tolerance: role of protein metabolism. Annu Rev Plant Physiol Plant Mol Biol. 1990;41:187–223.

    Article  CAS  Google Scholar 

  22. Thomashow MF. Plant cold acclimation: freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999;50:571–99.

    Article  CAS  PubMed  Google Scholar 

  23. Jaglo-Ottosen KR, Gilmour SJ, Zarka DG, Schabenberger O, Thomashow MF. Arabidopsis CBF1 overexpression induces COR genes and enhances freezing tolerance. Science. 1998;280:104–6.

    Article  CAS  PubMed  Google Scholar 

  24. Park SC, Lee CM, Doherty CJ, Gilmour SJ, Kim YS, Thomashow MF. Regulation of the Arabidopsis CBF regulon by a complex low-temperature regulatory network. Plant J. 2015;82:193–207.

    Article  CAS  PubMed  Google Scholar 

  25. Jia YX, Ding YL, Shi YT, Zhang XY, Gong ZZ, Yang SH. The cbfs triple mutants reveal the essential functions of CBFs in cold acclimation and allow the definition of CBF regulons in Arabidopsis. New Phytol. 2016;212:345–53.

    Article  CAS  PubMed  Google Scholar 

  26. Chinnusamy V, Ohta M, Kanrar S, Lee BH, Hong X, Agarwal M, et al. ICE1:a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. Genes Dev. 2003;17:1043–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Fursova OV, Pogorelko GV, Tarasov VA. Identification of ICE2, a gene involved in cold acclimation which determines freezing tolerance in Arabidopsis thaliana. Gene. 2009;429:98–103.

    Article  CAS  PubMed  Google Scholar 

  28. Doherty CJ, Van Buskirk HA, Myers SJ, Thomashow MF. Roles for Arabidopsis CAMTA transcription factors in cold-regulated gene expression and freezing tolerance. Plant Cell. 2009;21:972–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Chan ZL, Wang YP, Cao MJ, Gong YH, Mu ZX, Wang HQ, et al. RDM4 modulates cold stress resistance in Arabidopsis partially through the CBF-mediated pathway. New Phytol. 2016;209:1527–39.

    Article  CAS  PubMed  Google Scholar 

  30. Guo XY, Xu SJ, Chong K. Cold signal shuttles from membrane to nucleus. Mol Cell. 2017;66:7–8.

    Article  CAS  PubMed  Google Scholar 

  31. Kindgren P, Ard R, Ivanov M, Marquardt S. Transcriptional read-through of the long noncoding RNA SVALKA governs plant cold acclimation. Nat Commun. 2018;5:4561.

    Article  CAS  Google Scholar 

  32. Nichols PGH, Loi A, Nutt BJ, Evans PM, Craig AD, Pengelly BC, et al. New annual and short-lived perennial pasture legumes for Australian agriculture-15 years of revolution. Field Crops Res. 2007;104:10–23.

    Article  Google Scholar 

  33. Cook DR. Medicago truncatula - a model in the making! Curr Opin Plant Biol. 1999;2:301–4.

    Article  CAS  PubMed  Google Scholar 

  34. Tadege M, Ratet P, Mysore KS. Insertional mutagenesis: a Swiss Army knife for functional genomics of Medicago truncatula. Trends Plant Sci. 2005;10:229–35.

    Article  CAS  PubMed  Google Scholar 

  35. Young ND, Udvardi M. Translating Medicago truncatula genomics to crop legumes. Curr Opin Plant Biol. 2009;12:193–201.

    Article  CAS  PubMed  Google Scholar 

  36. Pennycooke JC, Cheng HM, Stockinger EJ. Comparative genomic sequence and expression analyses of Medicago truncatula and alfalfa subspecies falcate COLD-ACCLIMATION-SPECIFIC genes. Plant Physiol. 2008;146:1242–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Zhang LL, Zhao MG, Tian QY, Zhang WH. Comparative studies on tolerance of Medicago truncatula and Medicago falcata to freezing. Planta. 2011;234:445–57.

    Article  CAS  PubMed  Google Scholar 

  38. Zhang ZQ, Hu XN, Zhang YQ, Miao ZY, Xie C, Meng XZ, et al. Opposing control by transcription factors MYB61 and MYB3 increases freezing tolerance by relieving C-repeat binding factor suppression. Plant Physiol. 2016;172:1306–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Mohapatra SS, Wolfraim L, Poole RJ, Dhindsa RS. Molecular cloning and relationship to freezing tolerance of cold-acclimation-specific genes of alfalfa. Plant Physiol. 1989;89:375–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Wolfraim LA, Langis R, Tyson H, Dhindsa RS. cDNA sequence, expression, and transcript stability of a cold acclimation-specific gene, cas18, of alfalfa (Medicago falcata) cells. Plant Physiol. 1993;101:1275–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Zhao MG, Liu WJ, Xia XZ, Wang TZ, Zhang WH. Cold acclimation-induced freezing tolerance of Medicago truncatula seedlings is negatively regulated by ethylene. Physiol Plantarum. 2014;152:115–29.

    Article  CAS  Google Scholar 

  42. Tayeh N, Bahrman N, Sellier H, Bluteau A, Blassiau C, Fourment J, et al. A tandem array of CBF/DREB1 genes is located in a major freezing tolerance QTL region on Medicago truncatula chromosome 6. BMC Genomics. 2013;14:814.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Ewing B, Green P. Base-calling of automated sequencer traces using Phred. II Error probabilities. Genome Res. 1998;8:186–94.

    Article  CAS  PubMed  Google Scholar 

  44. Li S, Yu X, Lei N, Cheng Z, Zhao P, He Y, et al. Genome-wide identification and functional prediction of cold and/or drought-responsive lncRNAs in cassava. Sci Rep. 2017;7:45981.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Bigras FJ, Dumais D. Root-freezing damage in the containerized nursery: impact on plantation sites-a review. New For. 2005;30:167–84.

    Article  Google Scholar 

  46. Suzuki K, Nagasuga K, Okada M. The chilling injury induced by high root temperature in the leaves of rice seedlings. Plant Cell Physiol. 2008;49:433–42.

    Article  CAS  PubMed  Google Scholar 

  47. Hekneby M, Antolín MC. Sánchez-DíazM. Frost resistance and biochemical changes during cold acclimation in different annual legumes. Environ Exp Bot. 2006;55:305–14.

    Article  CAS  Google Scholar 

  48. Knight MR, Knight H. Low-temperature perception leading to gene expression and cold tolerance in higher plants. New Phytol. 2012;195:737–51.

    Article  CAS  PubMed  Google Scholar 

  49. Shu YJ, Liu Y, Li W, Song LL, Zhang J, Guo CH. Genome-wide investigation of microRNAs and their targets in response to freezing stress in Medicago sativa L., based on high-throughput sequencing. G3. 2016;6:755–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. McKersie BD, Chen Y, de Beus M, Bowley SR, Bowler C, Inze D, et al. Superoxide dismutase enhances tolerance of freezing stress in transgenic alfalfa (Medicago sativa L.). Plant Physiol. 1993;103:1155–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Sunkar R, Kapoor A, Zhu JK. Posttranscriptional induction of two cu/Zn superoxide dismutase genes in Arabidopsis is mediated by down regulation of miR398 and important for oxidative stress tolerance. Plant Cell. 2006;18:2051–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Swartzlander DB, McPherson AJ, Powers HR, Limpose KL, Kuiper EG, Degtyareva NP, et al. Identification of SUMO modification sites in the base excision repair protein, Ntg1. DNA Repair. 2016;48:51–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Chen HC, Chu P, Zhou YL, Li Y, Liu J, Ding Y, et al. Overexpression of AtOGG1, a DNA glycosylase/AP lyase, enhances seed longevity and abiotic stress tolerance in Arabidopsis. J Exp Bot. 2012;63:4107–21.

    Article  CAS  PubMed  Google Scholar 

  54. Jiang SY, Ramachandran S. Expansion mechanisms and evolutionary history on genes encoding DNA glycosylases and their involvement in stress and hormone signaling. Genome Biol Evol. 2016;8:1165–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Knight H, Mugford SG, Ulker B, Gao DH, Thorlby G, Knight MR. Identification of SFR6, a key component in cold acclimation acting post-translationally on CBF function. Plant J. 2009;58:97–108.

    Article  CAS  PubMed  Google Scholar 

  56. Avia K, Pilet-Nayel ML, Bahrman N, Baranger A, Delbreil B, Fontaine V, et al. Genetic variability and QTL mapping of freezing tolerance and related traits in Medicago truncatula. Theor Appl Genet. 2013;126:2353–66.

    Article  CAS  PubMed  Google Scholar 

  57. Guttman M, Rinn JL. Modular regulatory principles of large non-coding RNAs. Nature. 2012;482:339–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Gilmour SJ, Zarka DG, Stockinger EJ, Salazar MP, Houghton JM, Thomashow MF. Low temperature regulation of the Arabidopsis CBF family of AP2 transcriptional activators as an early step in cold-induced COR gene expression. Plant J. 1998;16:433–42.

    Article  CAS  PubMed  Google Scholar 

  59. Liu ZY, Jia YX, Ding YL, Shi YT, Li Z, Guo Y, et al. Plasma membrane CRPK1-mediated phosphorylation of 14-3-3 proteins induces their nuclear import to fine-tune CBF signaling during cold response. Mol Cell. 2017;66:117–28.

    Article  CAS  PubMed  Google Scholar 

  60. Quan MY, Chen JH, Zhang DQ. Exploring the secrets of long noncoding RNAs. Int J Mol Sci. 2015;16:5467–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Aird D, Ross MG, Chen WS, Danielsson M, Fennell T, Carsten R, et al. Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries. Genome Biol. 2011;12:1–14.

    Article  CAS  Google Scholar 

  62. Su ZQ, Labaj PP, Li S, Thierry-Mieg J, Thierry-Mieg D, Shi W, et al. A comprehensive assessment of RNA-seq accuracy, reproducibility and information content by the sequencing quality control consortium. Nature Biotechnol. 2014;32:903–17.

    Article  CAS  Google Scholar 

  63. Young ND, Cannon SB, Sato S, Kim D, Cook DR, Town CD, et al. Sequencing the gene spaces of Medicago truncatula and Lotus japonicas. Plant Physiol. 2005;137:1174–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357–U54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11:1650–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Sun L, Luo HT, Bu DC, Zhao GG, Yu KT, Zhang CH, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41:e166.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei LP, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35:W345–9.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40:D290–301.

    Article  CAS  PubMed  Google Scholar 

  70. Lin MF, Jungreis I, Kellis M. PhyloCSF: a comparative genomics method to distinguish protein coding and non-coding regions. Bioinformatics. 2011;27:i275–i82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics. 2008;9:559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  72. Pimentel H, Bray NL, Puente S, Melsted P, Pachter L. Differential analysis of RNA-seq incorporating quantification uncertainty. Nat Methods. 2017;14:687–90.

    Article  CAS  PubMed  Google Scholar 

  73. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11:R14.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  74. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank the technical support from the Test Centre of State Key Laboratory of Vegetation and Environmental Change.

Funding

This research was supported by National Natural Science Foundation of China (31671270) and Chinese Academy of Sciences (KFJ-STS-ZDTP-056).

Author information

Authors and Affiliations

  1. State Key Laboratory of Vegetation and Environmental Change, Institute of Botany, the Chinese Academy of Sciences, Beijing, 100093, People’s Republic of China

    Mingui Zhao, Tianzuo Wang, Tianyang Sun, Xiaoxi Yu, Rui Tian & Wen-Hao Zhang

  2. Research Network of Global Change Biology, Beijing Institutes of Life Science, the Chinese Academy of Sciences, Beijing, 100101, People’s Republic of China

    Mingui Zhao, Tianzuo Wang & Wen-Hao Zhang

  3. University of Chinese Academy of Sciences, Beijing, 100049, People’s Republic of China

    Tianyang Sun, Xiaoxi Yu, Rui Tian & Wen-Hao Zhang

  4. Inner Mongolia Research Center for Prataculture, the Chinese Academy of Sciences, Beijing, 100093, People’s Republic of China

    Wen-Hao Zhang

Authors
  1. Mingui Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Tianzuo Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Tianyang Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xiaoxi Yu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Rui Tian
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Wen-Hao Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

The authors have made the following declarations about their contributions: Conceived the research: Z. M. and Z. W-H. Designed the experiments: Z. M. Performed the experiments and analysed the data: Z. M., W. T., S. T., Y. X., and T. R. Wrote the paper: Z. M. and Z. W-H. All authors have read and approved the manuscript.

Corresponding author

Correspondence to Wen-Hao Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

Co-author Wen-Hao Zhang is an Associate Editor of this journal.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Statistical date of the RNA-Seq quality for experimental samples. Table S2. Correlation coefficient of lncRNAs for experimental samples. Table S3. Correlation coefficient of mRNAs for experimental samples. Table S4. Statistical date of the RNA-Seq results mapped to M. truncatula A17 genome for experimental samples. Table S5. GO enhancements of the putative targets of cold-responsive lncRNAs in leaves of M. truncatula seedlings. Table S6. GO enhancements of the putative targets of cold-responsive lncRNAs in roots of M. truncatula seedlings. Table S7. Primer sequences used for real-time quantitative RT-PCR.

Additional file 2: Fig. S1.

Density distribution of lncRNAs on eight chromosomes of M. truncatula seedlings.

Additional file 3: Fig. S2.

Number and density distribution of lncRNAs on chromosome of M. truncatula seedlings with and without cold treatment.

Additional file 4: Fig. S3.

Identification of mRNAs by high-throughput sequencing in M. truncatula seedlings.

Additional file 5: Fig. S4.

Characteristics of cold-responsive mRNAs in M. truncatula seedlings.

Additional file 6: Fig. S5.

Protein sequence alignment of the MtCBFs and AtCBFs was constructed by MAGE.

Additional file 7: Fig. S6.

Phylogenetic tree of MtCBFs and AtCBFs was constructed by MEGA.

Additional file 8: Data Set S1.

lncRNA TPM result.

Additional file 9: Data Set S2.

lncRNA sequence.

Additional file 10: Data Set S3.

mRNA TPM result.

Additional file 11: Data Set S4.

mRNA sequence.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhao, M., Wang, T., Sun, T. et al. Identification of tissue-specific and cold-responsive lncRNAs in Medicago truncatula by high-throughput RNA sequencing. BMC Plant Biol 20, 99 (2020). https://doi.org/10.1186/s12870-020-2301-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-020-2301-1

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us