Identification and functional prediction of cold-related long non-coding RNA (lncRNA) in grapevine

Plant long non-coding RNA (lncRNA) undergoes dynamic regulation and acts in developmental and stress regulation. In this study, we surveyed the expression dynamics of lncRNAs in grapevine (Vitis vinifera L.) under cold stress using high-throughput sequencing. Two-hundred and three known lncRNAs were significantly up-regulated and 144 known lncRNAs were significantly down-regulated in cold-treated grapevine. In addition, 2 088 novel lncRNA transcripts were identified in this study, with 284 novel lncRNAs significantly up-regulated and 182 novel lncRNAs significantly down-regulated in cold-treated grapevine. Two-hundred and forty-two differentially expressed grapevine lncRNAs were predicted to target 326 protein-coding genes in a cis-regulatory relationship. Many differentially expressed grapevine lncRNAs targeted stress response-related genes, such as CBF4 transcription factor genes, late embryogenesis abundant protein genes, peroxisome biogenesis protein genes, and WRKY transcription factor genes. Sixty-two differentially expressed grapevine lncRNAs were predicted to target 100 protein-coding genes in a trans-regulatory relationship. The expression of overall target genes in both cis and trans-regulatory relationships were positively related to the expression of lncRNAs in grapevines under cold stress. We identified 31 known lncRNAs as 34 grapevine micro RNA (miRNA) precursors and some miRNAs may be derived from multiple lncRNAs. We found 212 lncRNAs acting as targets of miRNAs in grapevines, involving 150 miRNAs; additionally, 120 grapevine genes were predicted as targets of grapevine miRNAs and lncRNAs. We found one gene cluster that was up-regulated and showed the same expression trend. In this cluster, many genes may be involved in abiotic stress response such as WRKY, Hsf, and NAC transcription factor genes.


Results
Data mining of transcriptome sequencing and identification of lncRNAs in grapevine. To systematically identify lncRNAs related to cold stress in grapevine, we performed whole transcriptome RNA-seq of grapevine cv. Cabernet Sauvignon that had been submitted to a cold-stress treatment of 4 °C. We generated an average of 12.65 gigabases (Gb) of raw reads per sample from the six samples used for Illumina RNA-sequencing. The total number of raw reads per control (CK) sample (plants were kept under a 16-h light/8-h dark photoperiod at 26 °C) ranged from 220842362 to 274931726, and the number of clean reads in each CK sample ranged from 216561108 to 270342092. The total number of raw reads in each cold treatment sample ranged from 191766324 to 233777742, and the number of clean reads in each cold treatment sample ranged from 186259776 to 223345340. The average mapping rate to the grapevine genome is 63.77%. In total, we identified 56732 transcripts, including 44644 known mRNA transcripts, 2031 known lncRNA transcripts, 7969 novel mRNA transcripts, and 2088 novel lncRNA transcripts. The transcripts of novel lncRNAs predicted here are listed in Table S1, and the transcript names and the related lncRNA gene IDs are listed in Table S2. The transcripts of novel mRNAs predicted here are listed in  Table S3, and the transcript names and the related mRNA gene IDs are listed in Table S4. In all samples, we identified 212 novel lincRNAs, 1 933 novel long intronic RNAs, and 688 novel lncNAT. We also found 1 893 known lincRNAs, 511 known long intronic RNAs, and 803 known lncNAT in total samples. In addition, we found that it was the most common for the lncRNAs to contain only one exon; lncRNAs containing two exons were the next most common, followed by lncRNAs containing three exons and four exons (Fig. 1A). The lncRNAs less than 500 bp long were most common, followed by the 500-1000 bp long lncRNAs and 1000-1500 bp long lncRNAs (Fig. 1B). We also found that most lncRNAs were located on chromosome 1 (Fig. 1C).
Variation in lncRNA expression among cold stress. In grapevine, 17 known lncRNAs were expressed only in the CK library and 97 known lncRNAs were expressed only in the cold-treated library. The expression heatmaps of all known and novel grape lncRNAs in the CK and cold treatment based on the Fragments Per Kilobase Million (FPKM) model are shown in Fig. 2A, and the box plot of expression levels of grape lncRNAs in the CK and cold treatment are shown in Fig. 2B. In both the control and cold treatments, the average expression level of the total lncRNAs was lower than that of the mRNAs in grapevine (Fig. 2C,D). Two-hundred and three known lncRNAs were significantly up-regulated (fold change > 2, P < 0.05) and 144 known lncRNAs were significantly down-regulated in cold-treated grapevine (fold change < −2, P < 0.05). In grapevine, VIT_203s0017n00360 was the lncRNA with the greatest increase of up-regulation by the cold treatment, followed by VIT_207s0031n00070 and VIT_201s0011n00530. VIT_209s0002n00340 was the lncRNA with the greatest down-regulation by cold treatment, followed by VIT_213s0158n00020 and VIT_213s0067n00110. These significantly up-and down-regulated lncRNAs were considered the differentially expressed known lncR-NAs (Fig. 3A, Table S5).
In grapevine, 17 novel lncRNAs were expressed only in the untreated library and 11 novel lncRNAs were expressed only in the cold-treated library. We identified 284 novel lncRNAs as significantly up-regulated (fold change > 2, P < 0.05) and 182 novel lncRNAs were significantly down-regulated (fold change > 2, P < 0.05) in cold-treated grapevine compared with in the CK. In grapevine, LXLOC_001173 was the lncRNA with the greatest up-regulation in the cold treatment compared with the CK, followed by LXLOC_004676 and LXLOC_028762. Compared with the CK, LXLOC_003867 was the lncRNA with the greatest down-regulation in the cold Prediction of target genes of cold-related lncRNA targets in cis-regulatory relationships. To investigate the possible functions of grape lncRNAs, we predicted the potential targets of lncRNAs in cis-regulatory relationships. We searched for known protein-coding genes located within 10 kb downstream and upstream of all the identified grape lncRNAs. These genes were thought to be the targets of lncRNAs in cis-regulatory relationships if the Pearson and Spearman correlation coefficients between the expression levels of these genes were ≥0.6 or ≤−0.6, and P < 0.05 41 .
Our results predicted a total of 2 527 target genes in cis-regulatory relationships of 1 650 lncRNAs in grapevine. In our study, significantly up-regulated or down-regulated lncRNAs were thought to be differentially expressed lncRNAs. Specifically, we found that 242 differentially expressed grapevine lncRNAs were predicted to target 326 protein-coding genes in cis-regulatory relationships, and many differentially expressed grapevine lncR-NAs targeted stress response-related genes such as CBF4 transcription factor genes, late embryogenesis abundant protein genes, peroxisome biogenesis protein genes, and WRKY transcription factor genes ( Table S6). The expression levels of some target genes in cis-regulatory relationships were positively related to lncRNAs. For example, VIT_216s0100n00030, LXLOC_027751, LXLOC_010422, and VIT_202s0025n00100 were up-regulated under cold stress compared to the CK. Based on our RNA-seq data, their target genes VIT_216s0100g00380 (CBF4 www.nature.com/scientificreports www.nature.com/scientificreports/ transcription factor), VIT_208s0058g00960 (transcription factor bHLH61), VIT_215s0046g02110 (late embryogenesis abundant protein Lea14-A), and VIT_202s0025g01280 (WRKY transcription factor 41) respectively, were also up-regulated under cold stress compared to the CK. The RNA-seq data was validated by the qRT-PCR results (Fig. 4). The expression levels of some target genes in cis-regulatory relationships were negatively related to lncRNAs. For example, compared to the CK, LXLOC_013001 was down-regulated under cold stress, but its target gene in the cis-regulatory relationship, VIT_217s0000g06350 (chlorophyll a-b binding protein 4), was up-regulated under cold stress when compared to the control. LXLOC_019156 was up-regulated under cold stress, but its target gene in the cis-regulatory relationship, VIT_202s0154g00610 (peroxisome biogenesis protein), was down-regulated under cold stress (Table S6, Fig. 5A). We calculated the correlation coefficient between the expression changes of lncRNAs and their target genes in cis-regulatory relationships under cold stress. As shown in Fig. 5B, the values of the x-axis are the log2fold change of lncRNAs (fold change = FPKM value of genes in the cold treatment/FPKM value of lncRNA genes in the control). The values along the y-axis are the log2fold change of their target genes in cis-regulatory relationships (fold change = FPKM value of genes in the cold treatment/FPKM value of genes in the control). The correlation coefficient was 0.53 (t-test, P < 0.05), indicating that the expression of overall target genes with a cis-regulatory relationship was positively related to the expression of related lncRNAs in grapevine under cold stress (Fig. 5B). The heatmap of expression of lncRNAs and their target genes in cis-regulatory relationships under cold stress based on the log2fold change value also showed that the expression of the overall target genes with cis-regulatory relationships were positively related to the expression of related lncRNAs in grapevine under cold stress (Fig. 5A).

Analysis of target genes of cold-related lncRNAs in trans-regulatory relationships.
To investigate the possible functions of grapevine lncRNAs, we predicted the potential targets of lncRNAs in trans-regulatory relationships. RNAplex software 42 was used to identify the lncRNA (parameters: >−30 binding energy) as described in a previous study 41 . The Pearson and Spearman correlation coefficients between the expression of these genes identified using RNAplex and the expression of related lncRNAs must be ≥0.6 or ≤0.6 and P < 0.05, or will be filtered out 41 .
In grapevine, we predicted a total of 574 target genes in trans-regulatory relationships with 422 lncR-NAs (Table S6). The results showed that 62 differentially expressed grapevine lncRNAs were predicted to target 100 protein-coding genes in trans-regulatory relationships such as NADH dehydrogenase subunit genes, UDP-glycosyltransferase genes, calcium-transporting ATPase genes, disease resistance protein genes, and glutamate receptor genes. However, most target genes in trans-regulatory relationships were unknown protein coding genes (Table S6). The expression levels of some target genes in trans-regulatory relationships were positively related with lncRNAs. For example, VIT_200s0225n00020 was down-regulated under cold stress, and its target gene, VIT_200s0246g00150 (NADH dehydrogenase subunit 5), in the trans-regulatory relationship was up-regulated under cold stress. Some target genes with trans-regulatory relationships were negatively related to lncRNAs (Fig. 5C).
We calculated the correlation coefficient between the expression changes of lncRNAs and their target genes in trans-regulatory relationships in the cold stress treatment. As shown in Fig. 5D, the values along the x-axis are the log2fold change of lncRNAs (fold change = FPKM value of genes in the cold treatment/FPKM value of lncRNA genes in the CK). The values along the y-axis are the log2fold change of their target genes in trans-regulatory relationships (fold change = FPKM value of genes in the cold treatment/FPKM value of genes in the control). The correlation coefficient was 0.71 (t-test, P < 0.05), indicating that the expression levels of overall target genes www.nature.com/scientificreports www.nature.com/scientificreports/ with trans-regulatory relationships were positively related to lncRNAs in grapevine under cold stress (Fig. 5D). The heatmap of expression of lncRNAs and their target genes in trans-regulatory relationships under cold stress based on the log2fold change value also showed that the expression of overall target genes with a trans-regulatory relationship were positively related to the expression of lncRNAs in grapevine under cold stress (Fig. 5C).
Validation of lncRNA expression using qRT-PCR. We performed qRT-PCR analyses to validate the RNA-seq results from six randomly selected grapevine lncRNAs, VIT_201s0010n00070, VIT_209s0002n00020, VIT_200s0179n00030, VIT_207s0141n00070, VIT_208s0007n00270, and VIT_207s0005n00480. The primers for qRT-PCR are listed in Table S8. The expression results were similar to the deep sequencing data (Fig. 7). VIT_200s0179n00030, VIT_207s0141n00070, and VIT_207s0005n00480 were shown to be up-regulated by the qRT-PCR data, showing a positive correlation with the deep sequencing results. VIT_201s0010n00070, VIT_208s0007n00270 and VIT_209s0002n00020 were down-regulated in both the qRT-PCR and RNA-seq results (Fig. 7, Table S5).
The relationships between grape mRNAs, lncRNAs, and miRNAs. We predicted the lncRNAs as targets or target mimics of miRNAs. In some previous studies, lncRNAs were found as both targets and target mimics of miRNAs 43,44 . As target mimics, lncRNAs could bind to miRNAs with a three-nucleotide bulge 43 . In our study, we only found lncRNAs that paired with miRNAs without any bulges. These lncRNA may be targets of miRNAs but not target mimics of miRNAs. Here, 212 lncRNAs as targets of miRNAs in grapevine were involved with 150 miRNAs (Table S9). Additionally, 120 predicted grapevine genes were both the target of grapevine miR-NAs and lncRNAs (Table S10).

Discussion
A previous study reported the existence of lncRNAs in plants 37 . As next generation sequencing technology developed, it became possible to identify lncRNAs including those identified in Arabidopsis, rice, maize, cassava 4,8,28,46 , and grapevine (http://genomes.cribi.unipd.it/DATA/V2/V2.1/lncRNA/). However, few studies have been conducted on the roles of lncRNAs involved in abiotic and biotic stress responses. In addition, there has been limited research conducted on the roles of lncRNA involved in abiotic stress response, such as response to cold stress, in grapevine. In this study, we detected the expression changes of lncRNAs in grapevine exposed to cold treatment www.nature.com/scientificreports www.nature.com/scientificreports/ and found 2 088 novel grapevine lncRNAs. Previous studies have also identified novel lncRNAs in other plant taxa including 6 500 novel lncRNAs in Arabidopsis thaliana 8 , 1 704 novel lncRNAs in maize 45 , and 682 novel lncRNAs in cassava 28 .
Here, we found that the average expression level of the total lncRNAs was lower than the average expression level of mRNAs in grapevine in both the control and cold treatment conditions (Fig. 2C,D). This indicates that the expression levels of total lncRNAs should be lower than mRNAs in grapevine. In A. thaliana, approximately 300 lncRNAs were evidenced to be differentially expressed under abiotic stressors 27,31 , and 318 cassava lncRNAs were differentially expressed under cold and drought conditions 28 . Here, we found 813 differentially expressed grapevine lncRNAs in the cold stress treatment, showing that more grapevine lncRNAs were differentially expressed under cold stress. We hypothesize that many grapevine lncRNAs may be related to cold stress and may play important roles in cold stress response. Though the expression levels in many lncRNAs changed in the cold treatment, the average expression levels of the total lncRNAs in the cold treatment were similar to the average expression levels of the total lncRNAs under control conditions (Fig. 2B).
We predicted the target genes of cold inducible grape lncRNAs, finding more target genes of cold inducible grapevine lncRNAs in cis-regulatory relationships than in trans-regulatory relationships. This indicated that the target genes in cis-regulatory relationships may be more related to cold stress response. We also analyzed the expression correlation between the total cold inducible grapevine lncRNAs and their target genes, and our results showed that the expression patterns were positively related.
The expression correlation between cold inducible grapevine lncRNAs and their target genes in trans-regulatory relationships were higher than in cis-regulatory relationships. However, some of the expression patterns of lncRNAs were negatively related to their target genes (Fig. 3). A previous study showed that lncRNAs could act as enhancers of gene expression 47 . In kiwifruit, the expression of both protein-coding genes and lncRNA  www.nature.com/scientificreports www.nature.com/scientificreports/ genes tended to be more positively correlated than negatively correlated in trans-regulatory relationships 48 . Here, we found that the expression of the overall target genes with a cis-regulatory relationship was also positively related to the expression of related lncRNAs in grapevine under cold stress.
Some target genes of cold inducible grapevine lncRNAs in cis-regulatory relationships may be involved in abiotic stress response such as VIT_216s0100g00380 (CBF4 transcription factor), VIT_215s0046g02110 (late embryogenesis abundant protein Lea14-A), and VIT_202s0025g01280 (WRKY transcription factor 41). These genes were also up-regulated in the cold stress treatment. Previous research has shown that CBF family genes   www.nature.com/scientificreports www.nature.com/scientificreports/ play critical roles related to control of an important pathway in the cold acclimation process 49,50 . CBF4 is one of the most important members for the over-wintering of grapevines 50 . Some LEA proteins have been shown to be involved in the freezing tolerance of plants 51 . Additionally, some WRKY transcription factors have been shown to be involved in modulating gene expression in plants during cold stress 52 . These cold stress-related genes were also up-regulated under cold stress. Therefore, we hypothesize that these cold stress-related genes could be regulated by related lncRNAs under cold stress. These lncRNAs may play important roles in cold stress tolerance and may be related to the regulation of these cold stress-related genes.
The GO analysis showed that the biological process terms that are related to cold stress lncRNAs contained the regulation of jasmonic acid (JA) mediated signaling pathway (GO: 2000022), regulation of defense response (GO: 0031347), regulation of signal transduction (GO: 0009966), hormone metabolic process (GO: 0042445), and regulation of hormone levels (Fig. 6A, Table S7). Jasmonic acid is related to cold stress response in plants 53 . Other hormones, such as abscisic acid (ABA), are related to abiotic responses in plants 54 . Molecular function terms of genes that are related to cold stress lncRNAs contained transcription factor activity, sequence-specific DNA binding (GO: 0003700), and transcription factor binding (GO: 0000989) (Fig. 6A, Table S7), indicating that many target genes were transcription factors or were related to transcription factors. These transcription factors may be involved in cold response and the regulation of other downstream genes involved in cold response. Cellular component terms of genes that were related to the cold-related lncRNAs contained photosystem (GO: 0009521) and photosystem I (GO: 0009535) (Fig. 6A, Table S7), showing that many target genes may be related to photosystems. Under cold stress, the photosystems have been shown to be related to cold tolerance 55 .
A previous study has shown that the lncRNAs that acted as target mimics could bind to miRNAs with three-nucleotide bulges 43 . However, our data did not predict similar target mimics that have been found in previous studies 43 , but the data did predict some targets that could bind to miRNAs without three-nucleotide bulges. These lncRNA may be targets of miRNAs. Similarly, a previous report has shown that lncRNAs acting as target genes could bind to miRNAs without bulges 44 . We found that lncRNAs and protein coding genes shared common miRNAs, which could target both lncRNAs and protein coding genes, and miRNAs and lncRNAs shared common target genes in grapevine. We hypothesize that the lncRNAs may regulate protein coding genes via complex pathways in grapevines.
The genes with the same expression trends were clustered together, and the genes in the same cluster may be involved in the same biological process 45 . We identified one cluster (Cluster 9) that showed an up-regulated expression pattern under cold treatment (Fig. 8). In this cluster, many genes may be involved in abiotic stress response such as WRKY transcription factor genes 51 , Hsf transcription factor genes 16 , and NAC transcription factor genes 56 . In Cluster 9, we also found 19 WRKY transcription factor genes, most of which were significantly up-regulated. Cluster 9 contained many lncRNAs and many protein coding genes that are the target genes of the lncRNAs in this cluster. Therefore, we suggest that the cluster may contain one or more pathways related to cold stress response and that lncRNAs may play important roles in cold stress response in this pathway. Because many WRKYs were found in Cluster 9, WRKY family members may play important roles in the key cold stress response pathway. Although none of the WRKY genes in Cluster 9 was a target gene of the lncRNAs, they may still be indirectly  www.nature.com/scientificreports www.nature.com/scientificreports/ regulated by lncRNAs or regulated by the expression of lncRNAs; however, this requires further study. Additionally, in this cluster, there are some genes related to anthocyanin or flavonoid biosynthesis such as VvMYBA1, flavanone 7-O-glucoside 2′-O-beta-L-rhamnosyltransferase, isoflavone-7-O-methyltransferase 9, and WD repeat-containing protein (Table 11) [57][58][59] . Previous studies have shown that abiotic stressors (such as cold or heat stress) may regulate anthocyanin or flavonoid biosynthesis-related genes 57,58 . Anthocyanins have been shown to be synthesized as protective compounds in response to cold stress 60 . Our cluster analysis showed that some key anthocyanin biosynthesis related genes may be located in pathways involved in cold stress response; therefore, lncRNAs in pathways involved in cold stress response are related to these anthocyanin biosynthesis related genes. Supporting our findings, previous studies have shown that biotic or abiotic stressors are related to the biosynthesis of anthocyanins or flavonols in grapevine 58,59,61 . Further studies should be conducted on the relationship between anthocyanin/ flavonoid biosynthesis pathways and cold stress as additional results will positively impact viticulture and breeding.

Materials and Methods
Plant materials. One-year-old self-rooted seedlings of the grapevine cv. Cabernet Sauvignon were grown and maintained in the greenhouse under a 16 h light/8 h dark photoperiod at 26 °C. For the cold stress treatment, plant materials under a 16-h light/8-h dark photoperiod were transferred to 4 °C for 4 hours. For the control (CK), plants were kept under a 16-h light/8-h dark photoperiod at 26 °C for 4 hours. The shoot apices with well-developed leaves from these plant materials were collected. Each treatment consisted of three independent replicates. RNA was isolated for the construction of RNA-seq libraries and real-time PCR analysis. T-test p-values < 0.05 are considered to be significantly changed, and "*" represents a p-value < 0.05.