Identification of long non-coding RNAs in response to downy mildew stress in grape

LncRNAs are transcripts with length greater than 200 nt and do not encode proteins, but they regulate gene expression. Transcripts over 200 nt in length and short ORFs lacking coding potential are considered as the core criteria for selecting lncRNAs from several transcription units (TU) of RNA-seq read. It has been reported that lncRNAs have many types according to different classification principles. According to the genomic location, it can be divided into three main types: intronic lncRNA, lincRNA, antisense lncRNA^[1], and also can be further subdivided into ten types^[2].

lncRNA also has many functions, which can be used as another classification standard, such as transcriptional interference, small RNA precursor, generate endo-siRNAs and so on^[3]. The expression level of lncRNA is much lower than that of protein-coding genes, but it has a high tissue specific or cell specific expression pattern^[4]. Therefore, it can also be classified and studied according to its expression specificity in various tissues^[2].

It has been widely reported that lncRNA participates in various life activities, but there are relatively few reports in plants. lncRNA play important roles in plant growth and development. Many lncRNAs in Arabidopsis have been proved to be involved in the growth and development process. COOLAIR and COLDAIR can be used as regulatory molecules of plant flowering inhibitor FLC^[5,6]. HID1 can function through PIF3, the key inhibitor of photomorphogenesis, and then inhibit the elongation and development of cR hypocotyl^[7]. In rice, long-day specific male sterility-related lncRNA (LDMAR) regulates photoperiod sensitive male sterility (PSMS)^[8]. In addition, lncRNA has been reported to play a certain role in cotton (Gossypium arboreaum) fiber development^[9], chickpea (C. arietinum) flower development^[10], mulberry (Morus notabilis) root^[11] and flower organ development and strawberry (Fragaria vesca) flower and fruit development^[12]. Recently, lncRNA has been studied more and more in plant disease resistance. Cui et al. found that lncRNA33732 activated by WRKY1 induced RBOH expression, which increased the accumulation of H₂O₂ in the early defense reaction of tomato against Phytophthora infestans^[13]. In addition, another lncRNA15492 was reported that its overexpression in tomato plants enhanced the resistance of tomato to Phytophthora infestans^[14]. In grapes, lncRNAs have been found to respond to Botrytis cinerea by being target genes or endogenous mimic targets of miRNAs^[15]. Meanwhile, the result of functional enrichment showed lncRNAs were related to PR-10 protein^[16]. It has been previously reported in Arabidopsis thaliana that lncRNAs can promote the expression of PR-1 protein to enhance plant disease resistance^[17]. Besides, other plant disease, such as Fusarium oxysporum of Arabidopsis^[18], cotton Verticillium wilt^[19], phytoplasma in paulownia tomentosa^[20], also involved the life activities of lncRNAs. However, the role of lncRNAs in downy mildew remains to be explored.

Grape downy mildew is a worldwide disease caused by Plasmopara viticola. It is widely distributed in all grape producing areas in China, and mainly harms the leaves. Different grape varieties had different resistance degrees to downy mildew. Some East Asian grapes show significant resistance to downy mildew, while Vitis vinifera shows susceptibility to downy mildew^[21]. Marutyan^[22] and Sriniasam & Jeyarajam^[23] found that the activity of polycasein oxidase was higher in grape downy mildew resistant hybrid varieties, but lower in susceptible varieties. The peroxidase activity in the affected area of grape leaves increased and necrotic spots were formed^[24]. Wu et al.^[25] obtained some genes and pathways that might be related to downy mildew resistance by Solexa sequencing. The enrichment of differentially expressed genes (DEGs) involved many metabolic process, such as ribosomal, photosynthesis, amino acid and sugar. Genetic and expression analysis showed that PR protein coding genes^[26], NPR1 homologous genes^[27], Rpv1^[28] and glyoxalase-like gene families^[29] were related to the resistance mechanism of downy mildew in grapevines. RNAseq analysis of V. pseudoreticulata after infection of P. viticola revealed a molecular mechanism that oxidative stress and cell wall biosynthesis and modification were important to resist P. viticola^[30]. Also, some genes like VpPR10.1^[31], VvPR4b^[32], VvPYL4^[33] had been shown to play a role in resistance to grape downy mildew.

The molecular mechanism of the interaction between Plasmopara viticola and grapevine has been widely reported. During infection, oomycetes secrete effector proteins that typically suppress plant immunity by triggering sensitivity (ETS)^[34]. The two main types of secretory effectors identified in the oomycetes' genome are RxLR and CRN^[35]. Certain oomycetes effectors with specific protein motif RxLR are recognized by plant resistance genes, and this interaction triggers immunity (ETI), leading to local cell death or hypersensitivity (HR)^[36], they even inhibit plant RNA silencing^[37]. Breakthroughs have been made in the research on the mechanism of grape downy mildew resistance and resistance gene proteins, but the relevant regulation mechanism of grape lncRNA against downy mildew resistance is still unclear. To better understand the network of grape resistance, an investigation of dynamic changes in gene expression including lncRNAs during the course of disease is needed.

In this study, we identified 4,011 lncRNAs that were expressed in grape leaves of two species (Vitis piasezkii accession Liuba-8 (LB), V. vinifera cultivar Pinot Noir (PN)) at six stages after Plasmopara viticola inoculation by using RNA-seq. Next, we analyzed the characteristics of grape lncRNAs, including chromosome distribution, length, number of exons, classification and expression pattern. Moreover, we predicted cis-trans target protein-coding genes for lncRNAs and explored their possibility to be endogenous pseudo target mimics (eTMs) of known miRNAs. Finally, a lncRNA which may play a positive role in the early stage of grape downy mildew resistance was verified by transient transformation. Overall, our study provided a preliminary understanding of lncRNAs of grape, and helped to assess the functions of lncRNAs in plant disease resistance.

Grape plants 'Liuba-8' and 'Pinot Noir' are all from the Grape Planting Resource Garden of College of Horticulture, Northwest A&F University (Yangling, China). In the season from June to September, leaves with oil spot symptoms were collected from grape resource nursery, sporangium was collected by indoor culture, and sporangium suspension concentration was adjusted to 5 × 10⁴ mL⁻¹ to inoculate isolated leaves. The fully expanded leaves third to the sixth from the top were taken back to the laboratory in an ice box. The leaves were collected at 0, 12, 24, 48, 96, and 120 hpi after inoculation, with 0 hpi as the control samples. These materials were frozen in liquid nitrogen immediately after harvest and stored at −80 °C for subsequent analysis.

Transcriptome sequencing and identification of lncRNAs

Total RNA of grape leaves was extracted according to the instructions of the Plant RNA Kit (OMEGA), and the concentration and purity of RNA were determined by spectrophotometer. The extracted total RNA was sent to Beijing novogene for cDNA library construction and sequencing. Sequencing libraries were produced using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) and index codes were added to attribute sequences for each sample. The clustering of the index-coded samples was carried out on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina). Then, the library preparations were sequenced on an Illumina Hiseq platform, and paired end reads were generated. Hisat2 compared transcript reads of transcriptome sequencing results to genome, and assembled transcripts with Stringtie to obtain 58,489 new transcripts for lncRNAs screening; Genes that were longer than 200 bp were assembled; Compared with NR protein database, homologous search filtering was carried out, and the sequence without coding conservative protein domain was retained; Coding Potential Calculator (CPC) software was used to evaluate the protein coding potential of genes^[38,39] to screen possible non-coding RNAs; Housekeeping RNAs (e.g. rRNA, snRNA, tRNA, snoRNA) was filtered from Rfam database. Finally, we found highly reliable hypothesized lncRNAs in grapes.

Localization of lncRNAs and protein coding genes in grape genome

Circos(circos-0.69-9)^[40] was used to map the location of lncRNAs and protein coding genes in grape, and the relevant information of grape chromosome came from NCBI grape genome (www.ncbi.nlm.nih.gov/genome/401).

Classification of lncRNAs

According to the relative position with their nearest protein coding genes, the annotated lncRNAs are divided into four categories: (i) Intergenic lncRNAs (lincRNAs): lncRNAs do not have any overlap with other protein coding genes; (ii) Intron lncRNAs: lncRNAs were transcribed from the intron part of coding genes; (iii) sense lncRNAs: lncRNAs were overlapping with exons of protein-encoded transcripts on the same chain; (iv) Antisense lncRNAs: lncRNAs were overlapping with exons of protein-encoded transcripts on another chain. Perl scripts were developed to classify these four categories and is available on Github (https://github.com/Zhiqiang-hao/Script-for-lncRNAs-classification).

Transcription length and exon number distribution of lncRNAs and protein coding genes in grape

LncRNAs and protein coding genes were analyzed from the perspective of transcription length and exon number. Transcription lengths were classified as ≤ 300, 300−400, 400−500, 500−600, 600−700, 700−800, 800−900, 900−1,000 and > 1,000 nucleotides. Exon numbers were classified as 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 and > 10. Then, we calculated the proportion of different types of lncRNAs and protein coding transcripts using Excel.

Cis-trans regulating target gene prediction

For the prediction of Cis target gene, it is considered that the function of lncRNA is related to the protein coding gene near its coordinates. At first, the gene encoding protein located near its upstream or downstream (10 kb) was found, and then the target gene was enriched by GO function to preliminarily predict the possible function of lncRNAs. GO annotations and enrichment were using R packages. Using LncTar package^[41] and WGCNA to predict the trans-regulated coding genes

Prediction of miRNA target genes and endogenous target mimics from lncRNAs

Using perl language, according to the methods used in previous research^[42], grape lncRNA, which is the target gene of miRNA, is predicted as follows: (i) At most, only one mismatch is allowed from 9 to 12 positions from the 5' end of miRNA; (ii) The protrusion or mismatch in other areas cannot exceed 4 nt; (iii) Continuous mismatch is not allowed. The perl script is available on Github (https://github.com/Zhiqiang-hao/Script-for-miRNA-target-predict).

According to Wu et al.^[43], the prediction method of plant lncRNAs as the endogenous target genes of microRNA was predicted. Studies have shown that the double-stranded RNA formed by RNA as an endogenous pseudo-target gene and its miRNA often has a 3-nucleotide protrusion in the middle of the binding site, which is considered to be related to the loss of the cleavage function of miRNAs^[44,45]. In this prediction method, we require that: (1) The protrusion of three nucleotides of endogenous pseudo target gene should be at the 9−12 positions starting from the 5' end of miRNA sequence; (2) The process of endogenous pseudo target gene consists of only three nucleotides; (3) The 2−8 positions from the 5' end of the miRNA sequence should be completely complementary and paired; (4) There should be no more than three mismatches and G/U pairs between endogenous pseudo target genes and miRNA, except for core processes.

The Arabidopsis eTMs sequences of miR156, miR160, miR171 and miR172 were downloaded from Wu et al.^[43]. Clustalw (www.genome.jp/tools-bin/clustalw) and weblogo (http://weblogo.threeplusone.com/) are used to generate multiple sequence comparison and sequence identification of eTMs of grape and Arabidopsis.

Analysis of lncRNA expression

The FPKM (Fragments Per Kilobase Million) value is used to reflect the expression abundance of Unigenes. According to the calculation standard of general transcriptome sequencing, FPKM value ≥ 0.1 indicates that genes are expressed. Venn diagrams used local script R language and online website drawing (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Selecting FDR ≤ 0.05 and P-value < 0.01 as differentially expressed genes, FC < −1 as down-regulated expression and FC > 1 as up-regulated expression (FC is log fold change to the base 2). The differentially expressed genes at 0.5, 1, 2, 4 and 5 dpi compared with 0 h were selected as the differentially expressed genes possibly related to grape downy mildew.

Using limma packet in R packet to calculate Pearson coefficient of MSTRG.12742.1 and protein-coding genes, considered |PCC| > 0.8 and P-value < 0.001 as co-expressed genes by trans-regulation.

PN leaves transient experiment

RNA was extracted and reversely transcribed into cDNA with random primers. The gene was amplified from LB cDNA with PCR, and then connected to the pCambia2300-GFP vector which was modified previously by our laboratory personnel. Agrobacterium strain GV3101 containing p2300-GFP-lncRNA vectors were grown at 28 °C. After shaking, cultures were harvested and re-suspended in infiltration buffer to a final OD₆₀₀ of 0.6 and left at room temperature for 3 h. The isolated leaves were treated with sterile needles and then vacuumed in the Agrobacterium suspension. The PN leaves infiltrated with empty vector was used as control. Each inoculation was carried out three times, and each time more than three different plants were infiltrated. The samples were incubated for 24 h in darkness at 22 ± 2 °C with 100% RH, sporangial suspension was inoculated and followed by incubation at the same temperature under a 16-h light/8-h photoperiod^[46,47]. Samples were taken at 1, 3, 5 and 7 d after inoculation stained with aniline blue and observed under epifluorescence UV microscopy.

The plants we used for RNA-seq are 'Liuba-8'(LB) and 'Pinot noir'(PN), which are resistant and susceptible to downy mildew, respectively. The transcription reads of the RNA-seq results were compared to the genome by Hisat2, and 58,489 new transcripts were assembled by Stringtie. In order to identify candidate lncRNAs, 58,489 transcripts were screened and filtered according to five strict criteria. First, the sequence length was ≥ 200 nt, and 58,487 transcripts were filtered. Long non-coding RNA was generally believed to have short ORF but not able to encode polypeptides longer than 100 amino acids^[48−51], and 13,152 transcripts with ORF < 300 were selected. By comparing with the NR protein database, homology search and filtration were carried out on the remaining 13,152 transcripts, and 4,325 sequences that did not encode conserved protein domains were retained. Next, protein coding potential was assessed using CPC to eliminate 260 possible coding transcripts. Besides, housekeeping ncRNA (rRNA, snRNA, tRNA, snoRNA) was filtered out in the Rfam database, and 4,011 grape lncRNAs were obtained (Fig. 1).

Figure 1. 

Process chart of screening lncRNAs candidates from transcriptome sequencing. The transcription reads of the RNA-seq results were compared to the genome by Hisat2, and new transcripts were assembled by Stringtie. Transcripts were filtered with the five criteria for identification of lncRNAs. (1) The sequence length was ≥ 200 nt; (2) ORF < 300; (3) Non encoded conserved protein domains; (4) Little protein coding potential; (5) Not housekeeping RNA.

Characterizations of lncRNAs combine with protein-coding genes in grapevine

We mapped the lncRNAs and annotated genes onto V. vinifera chromosomes, and the obtained Circos diagram clearly showed that lncRNAs in grapes were not uniformly distributed on chromosomes (Fig. 2a). Similar to the protein-coding genes, the gene density near telomeres was significantly higher than the heterochromatin area around centromeres. This result suggested that lncRNAs may have similar transcription characteristics with protein coding genes in grape.

Figure 2. 

Characterizations of lncRNAs and protein-coding genes in Liuba-8 and Pinot Noir. (a) Genome-wide distribution of grape lncRNAs compared with that of protein-coding genes (RNA-seq). Chromosomes are shown with different colours and in a circular form as the outer thick track. The inner chromosome scale (Mb) is labelled on each chromosome. For the inner two tracks, the abundance of lncRNAs and protein-coding genes in physical bins of 0.5 Mb for each chromosome are shown by blue and grey columns, respectively. (b) The number of lncRNAs on each chromosome. (c) Classification of grape lncRNAs according to its genomic position and overlap with protein-coding genes.

According to the relative position of the nearest protein coding gene on the chromosome, we divided lncRNAs into four types: sense lncRNA, antisense lncRNA, intron lncRNA and intergenic lncRNA. There were also a few lncRNAs of unknown classification, which have a large span and span at least two classification ranges, including some regions in the interzone and some in intron or antisense intron. It can be seen from Fig. 2b that the majority of lncRNAs belong to the intergenic, which is consistent with previous studies^[40]. In addition, a lot of lncRNAs were antisense genes of protein coding genes. Consistent with previous studies, there were many antisense lncRNAs in Arabidopsis thaliana^[52]. At the same time, the numbers of four types of lncRNAs from Watson (−) chain and Crick (+) chain were similar (Fig. 2c), which is consistent with the previous research results in tomato^[53].

In previous studies, lncRNAs are shorter than protein-coding genes and have fewer exons in plants and animals^{[50,51,53,54]}. To determine whether grape lncRNAs have these characteristics, we analyzed the length and exon number distribution of 4,011 lncRNAs, and compared them with all grape predictive protein coding genes (26,346). These 4,011 lncRNAs ranged from 200 bp to 7,047 bp in length, with an average length of 582 bp. We found that about 88% of lncRNAs were between 200 bp and 1,000 bp in size, and only 12% were larger than 1,000 bp. In contrast, about 85% of transcripts of protein coding genes were larger than 1,000 nucleotides (Fig. 3a). The results showed that most (90%) genes encoding grape lncRNAs only contained one or two exons, while the number of exons of genes encoding proteins ranged from 1 to ≥ 10 (Fig. 3b). These results indicated that unlike protein coding genes, most grape lncRNAs were relatively shorter and contain fewer exons.

Figure 3. 

The length and exons of lncRNAs and protein-coding transcripts. The distribution of (a) length and (b) numbers of exons of lncRNAs in comparison with protein-coding transcripts of grape (Vvinifera_145_Genoscope.12X.gene).

Prediction of cis-regulated target protein-coding genes of lncRNAs

The function of lncRNAs is related to the protein coding genes near its coordinates, and the lncRNAs located upstream and downstream of the coding protein may intersect with the promoter or other cis-acting elements of the coexpression genes, thus regulating the gene expression at the transcription or post-transcription level^[55].

The gene encoding the proximal protein is located within 10 kb, which was searched as the target gene of cis-regulation of lncRNAs. A total of 3,643 lncRNAs were found to have potential cis-regulatory effects on 6,622 protein coding genes. Among these, 70% lncRNAs have more than one target gene, of which 60% have two to four target genes, and only one lncRNAs has 10 target genes (Fig. 4a). Up to 74% of protein coding genes only correspond to one lncRNA, and only two protein coding genes are cis-regulated by up to 17 lncRNAs (Fig. 4b).

Figure 4. 

LncRNAs and their cis-regulated target protein-coding genes in Liuba-8 and Pinot Noir. (a) The number of target protein-coding genes cis-regulated by lncRNAs. (b) The number of lncRNAs that have potential cis-regulatory effects on protein-coding genes

Next, to have a preliminary understanding of the main functions of lncRNAs, the GO functional enrichment of the cis-regulated coding target genes of DElncRNAs was carried out. The results showed that these genes were involved in many biological processes in vivo (Supplemental Fig. S1). Many galactosyltransferase and inhibitor activities are the more representative groups in the molecular functions. In the category of cellular components, the transcripts corresponding to histone acetyltransferase (HATs) complexes are the most abundant. Besides, metabolic processes like galactose metabolic process and hexose metabolic process are the most abundant groups in the category of biological process, indicating that extensive metabolic activity occurred in plants treated with downy mildew.

Prediction of trans-regulated target protein-coding genes of differentially expressed lncRNAs

lncRNA can also regulate mRNA expression through trans-regulation. The function of lncRNA is related to their co-expressed protein-coding genes rather than the position of coding genes. The prediction of trans regulation depends on the interaction between lncRNAs and related mRNAs through the formation of complementary hybridization. Then through the co-expression association to screen the genes with relevant expressions and predict them as possible target genes, forming a co-expression network with several core genes as the center. The binding energy of mRNA and 91 DElncRNAs were calculated by LncTar software^[41]. The WGCNA (Weighted Gene Co-Expression Network Analysis) of the expression level of these genes was further carried out, and the top 25% of variance was screened, which was divided into 18 modules, of which seven modules had no lncRNA (Fig. 5a). Taking lncRNAs with high expression as the network center, some of the modules was constructed (Fig. 5b). The GO enrichment of these modules showed the most abundant group in biological processes is metabolic process, which produced many metabolites to improve plant disease resistance (Supplemental Fig. S2). The module with MSTRG.12742.1 that highly expressed in LB not PN after infection was shown mainly enriched in iron-sulfur cluster binding which has been demonstrated to be involved in grape immunity to downy mildew^[56] (Supplemental Fig. S3).

Figure 5. 

WGCNA of differentially expressed lncRNAs and its trans target genes. (a) Cluster Dendrogram of the DElncRNA and mRNA. (b) The network construction of lncRNAs and its trans-target genes.

Prediction of lncRNAs as target genes and endogenous target mimetics of miRNAs

The interaction between miRNAs and lncRNAs is the important way for lncRNAs to play its regulatory role^[57]. The lncRNAs can be targeted by miRNAs^[51], and can also act as eTMs (endogenous mimic targets) of miRNAs^[49]. To check whether lncRNA is the real target of miRNAs, we examined 4,011 lncRNAs and 184 V. vinifera miRNAs using local scripts, and found that 117 grape miRNAs targeted 184 lncRNAs (Supplemental Table S1). Meanwhile, six lncRNAs were predicted to be the potential eTMs for 15 miRNAs (Table 1 and Supplemental Table S2). It can be found that the interaction between lncRNAs and miRNAs may be widespread in grape. Therefore, the research progress of miRNAs related to grape have important reference value for the study of grape lncRNAs.

Table 1.  lncRNAs are acting as endogenous target mimics (eTMs) of miRNAs (miRNA156, miRNA160) in Vitis vinifera.

lncRNAs	Sequences	miRNA
MSTRG.14598.2		vvi-miR156e
MSTRG.13285.1		vvi-miR160c
MSTRG.13285.1		vvi-miR160d
MSTRG.13285.1		vvi-miR160e

In a previous study, several eTMs of A. thaliana miRNAs were screened, including miRNA156, miRNA159, miRNA160, miRNA166, miRNA167, miRNA171, miRNA172, miRNA395 and miRNA399^[43]. Sequence comparison between the selected eTMs and those reported in Arabidopsis revealed that the eTMs of the same miRNAs were conserved in both species (Fig. 6). The nucleotides that form the bulge region are not conserved in different species, or even in the same species.

Figure 6. 

Syntenic analysis of eTM156s, eTM160s,and eTM166s between V. vinifera and A. thaliana. (a) Sequence alignment of eTMs for miR156 (eTM156s). The bases with a colored background are the sequences paired to miR156. (b) Sequence alignment of eTMs for miR160 (eTM160s). The bases with a colored background are the sequences paired to miR160. (c) Sequence alignment of eTMs for miR166 (eTM166s). The bases with a colored background are the sequences paired to miR166.

Expression patterns of lncRNAs in grape leaves after inoculation with grape downy mildew

According to the general screening criteria for RNA-seq expression, among the 4,011 lncRNAs screened, the FPKM value in each sample was between 0 and 0.1, which was deemed not to be expressed. Among the expressed lncRNAs, the expression level of about 2/3 lncRNAs is extremely low in that the FPKM value is concentrated between 0.1 and 1 (Table 2). The largest FPKM value of LB and PN are 608.978 and 877.953 respectively. By analyzing the expression of lncRNAs in these two varieties, it was found that the expression number of lncRNAs in each period was above 2,000, and the most lncRNAs expressed in each variety were LB1 (2,542) and PN5 (2,190). 3.81% (153/4,011) lncRNA was only expressed in a specific period of a certain variety, for example, some lncRNAs appeared only in LB leaves inoculated with downy mildew for 1 d. Among them, at 0 dpi after inoculation of LB leaves, the most lncRNAs were specifically expressed, up to 27 (Table 3).

A total of 3,638 lncRNAs were expressed in grape leaves within 5 d after downy mildew inoculation, and 2,263 lncRNAs were expressed in both varieties (Fig. 7a). We identified 774 and 601 lncRNAs specifically expressed in LB and PN leaves after inoculation, respectively. Most lncRNAs are selectively expressed in some stages, and a considerable number of lncRNAs are only expressed in one stage. The largest amount of co-expression of lncRNA was at 0 dpi, and the least was at 2 dpi (Fig. 7b). The LB-specific expression of lncRNAs increased at 1 dpi, while PN decreased, with the largest difference at 2 dpi. At each stage, the expression number of LB total lncRNAs and specifically expressed lncRNAs was always higher than the number of lncRNAs expressed by PN. The overall trend in the amount of lncRNA expression was consistent with the progress of P. viticola infection. All of the above results indicate that lncRNAs probably play a role in grape disease resistance. In total, 3,037 lncRNAs were expressed within the five days of inoculation in LB leaves (Fig. 7a), among which 1856 lncRNAs were expressed throughout the six stages (Fig. 7c). In the PN, a total of 2,864 lncRNAs were expressed in leaves, of which 1,447 were continuously expressed during the entire developmental process (Fig. 7d).

Figure 7. 

The expression of lncRNAs in grape. (a) Co-expression of lncRNAs in leaves of Liuba-8 and Pinot Noir after inoculation with downy mildew. (b) The expression difference of lncRNAs in leaves of Liuba-8 and Pinot Noir inoculated with downy mildew at the same period. (c) A summary of the number of lncRNAs expressed in different stages of Liuba-8 inoculated with downy mildew. (d) A summary of the number of lncRNAs expressed in different stages of Pinot Noir inoculated with downy mildew.

The function analysis of MSTRG.12742.1

According to the screening criteria of differentially expressed genes, lncRNAs significantly up-regulated and down-regulated in two varieties were screened respectively (compared with 0 d). The result showed that the up- and down-regulated genes in LB were both more than genes in PN (Supplemental Fig. S4). The genes with high FPKM, which means high expression than genes in PN (Supplemental Fig. S4). The genes with high FPKM, which means high expression level, were selected from differentially expressed genes for qRT-PCR verification (Supplemental Fig. S5). The results showed that the expression level of MSTRG.12742.1 increased significantly in LB and much higher than PN (Fig. 8), which was consistent with transcriptome data. It was predicted that MSTRG.12742.1 might play a positive role in grape downy mildew resistance.

Figure 8. 

Expression pattern of MSTRG.12742.1 in V. piasezkii accession Liuba-8 and V. vinifera cv. Pinot Noir. Expression pattern of MSTRG.12742.1 at 0, 6, 12, 24, 48, 72, 96 and 120 hpi in (a) V. piasezkii accession Liuba-8 and (b) V. vinifera cv. Pinot Noir.

Pearson correlation coefficient is used to calculate the correlation between MSTRG.12742.1 and the amount of protein coding gene expression, the genes which | had a Pearson | > 0.4 and P-value < 0.01, were selected as possible target genes. The highest positive correlation was VIT_204s0008g02671.1 (COR = 0.943) (Fig. 9a). GO function enrichment analysis was carried out for the predicted target genes, and the bubble diagram showed the top 15 items of enrichment in each classification. It was found that the target genes were mainly enriched in the photosynthetic membrane (Supplemental Fig. S6), which may be involved in the related process of light reaction.

Figure 9. 

The verification of MSTRG.12742.1 and its predicted target gene VIT_204s0008g02671.1. (a) The Pearson correlation coefficient of the expression of MSTRG.12742.1 and VIT_204s0008g02671.1. (b) MSTRG.12742.1 and predicted target genes expression after instantaneous transformation of Vitis vinifera cv. Pinot Noir leaves.

The expression level of MSTRG.12742.1 was significantly increased in the leaves of the transiently recombinant vector (relative to the leaves of the wild-type and transformed empty vector) (Fig. 9b). According to Pearson correlation coefficient, the expression of the protein-coding gene with the highest correlation was examined. The result showed that there was significant increase in VIT_204s0008g02671.1, which was same as MSTRG.12742.1.

Identification of lncRNAs in response to downy mildew stress

To further determine whether MSTRG.12742.1 has a positive effect on grape downy mildew resistance, the leaves over-expressed genes were inoculated with P. viticola, and the P. viticola colonization was examined (Fig. 9b). Under epifluorescence UV microscopy, it was observed that the leaves over-expressed genes had only substomatal vesicles at 1 dpi, while the leaves over-expressed empty vector had already appeared primary hyphae (Fig. 10a, e); At 3 dpi, a small amount of hyphae extended in the tissue, and the hyphae grew relatively fast in no-load (Fig. 10b, f); Mycelial network was formed at 5 dpi (Fig. 10c, g), and a large number of sporangiophore extended from stomata at 7 dpi and sporangia was distributed at the top (Fig. 10d, h). By comparison, it is found that there are obvious differences in P. viticola colonization only in the early stage, and the differences become smaller in the later stage, even though there is almost no obvious difference. It is preliminarily predicted that this gene may play a role in the initial stage of infection by P. viticola.

Figure 10. 

P. viticola development in leaf discs of transient expression of MSTRG.12742.1 and empty vector as revealed by staining with aniline blue. (a1)−(h1) is stained with aniline blue under epifluorenscense UV-microscopy. (a2)−(h2) is under bright-field illumination in the field of vision corresponding to (a1)−(h1). St stomata, sv substomatal vesicle, ph primary hyphae, hy hyphae, sph sporangiophores, sp sporangia

Non-coding RNA was first considered as the noise of genome transcription. Recently, many studies have confirmed that lncRNA plays important roles in many life activities and has become one of the most popular research fields. With the development of technologies, the mechanism of action of lncRNA in animals has been widely studied, while the research on plant lncRNA is still in its infancy, even in model plants. Grapevine is one of the economically important fruit crops in the world, and there are many studies on its disease resistance, but there are few reports on lncRNA. Here, on the basis of transcriptome sequencing, we systematically analyzed the expression of lncRNAs in grape leaves infected by P. viticola at different stages, explored the characteristics of lncRNAs in grape plants, and compared with coding genes, preliminarily identified the lncRNAs related to grape downy mildew and its expression characteristics.

Genome-wide transcriptome analysis combined with computational methods is widely used to discover new lncRNAs in animals and plants. Because lncRNA is poorly conserved in plant species including model plants such as Arabidopsis thaliana and rice, it is impossible to use these resources to identify homologous lncRNA in other plant species^[58,59]. By using RNA-seq and new screening strategies to re-identify lncRNA in grape will offer new annotation information to the grape genome.

Comparing with annotated genes, in previous studies, lncRNAs possess short length and fewer exons in both plants and animals^{[50，51,53,54]}. In the grapevine, this characteristic also exists, which may reveal that this is an inherent characteristic of lncRNAs regardless of plant species and may contribute to its non-ability of coding protein. Moreover, the location of lncRNAs is related to its function and most grape lncRNAs belong to long intergenic non-coding RNA, consistent with other species. Because it is located between coding genes, it is easier to study and obtain than other types and there are also many reports. Above all, it is probably the best breakthrough point for grape lncRNAs research.

The target genes of grape lncRNAs regulate the related activities of plant immunity

lncRNA do not have the ability to encode proteins and cannot directly translate proteins to regulate plant growth and development. Therefore, it generally plays its role by regulating other genes. lncRNA could regulate gene expression through cis-acting mechanisms. cis-regulation refers to the regulation of neighboring protein coding genes located on the same chromosome^[60]. We searched the genes located within 10 kb of lncRNAs as the target gene by cis-regulation. The results showed most grape lncRNAs had target coding genes, even more than one (Fig. 4). lncRNA could cis-regulate the transcription of adjacent protein coding genes, precisely control the spatial and temporal expression of these developmental diversity gene sites, and participate in related development and other biological processes^[61]. The main functions of lncRNAs can be predicted by analyzing the functional enrichment of the target genes, which lays a foundation for subsequent cis analysis. The GO enrichment results showed that they were mainly enriched in the groups about glycosyltransferase (GT) and HATs, which are crucial for plant immunity. As reported, GT can respond to plant disease resistance by modifying salicylic acid and jasmonic acid^[62]. HATs plays a significant role in plant immunity, which can induce the expression of resistance genes, such as PR1, PR5 and WRKY18^[63]. To verify whether the lncRNAs could contribute to grape to resistant biological stress like downy mildew, further study on cloning and functional annotation of lncRNAs are needed. The trans-regulated target genes of lncRNAs refer to genes that are far away from them in position and regulate mRNA expression through sequence complementary pairing, and the expression level of the two genes will show a certain correlation. In this study, the binding energy of lncRNAs and mRNAs was calculated by Perl language and lncTar package^[41]. Then, WGCNA is used to construct the network of gene expression related modules^[64]. Due to the large amount of data to calculate the binding energy of all genes, we selected DElncRNAs that are more likely to be associated with grape downy mildew for prediction. As shown, the function of those trans-target genes were mainly enriched in metabolic processes which is associated with plant disease resistance.

lncRNA works with miRNA to regulate plant disease resistance

miRNA participates in the regulation of plants in the process of plant disease resistance. The research in transgenic tobacco showed that overexpression of miR396a-5p in tobacco could significantly reduce the resistance to Phytophthora nicotianae^[65]. miRNA is involved in grape response to stress. For example, the expression of miR156, miR162, miR166 and miR167a is related to grape leaf curl disease^[66]. In addition, some miRNAs were screened out and might be involved in the process of resistance to white rot of V. davidii^[67]. The lncRNAs can also function through miRNA. It was found that when lncRNA23468 was overexpressed in tomato, the expression of miR482b and its target gene NBS-LRR decreased and increased respectively, and the disease resistance was enhanced^[68]. lncRNAs play a role as the target genes of miRNAs. lncRNA has introns and other fragments, which reach thousands of bases in length, providing a good material basis for adsorption and binding of a large number of miRNAs. By competing to occupy a large number of miRNAs in the cell, it reduced its ability to interfere with the protein encoded by the target gene. Many kinds of miRNAs integrated regulation modes of lincRNA were detected in humans, some of which are important in breast cancer^[69]. The lncRNA would be as an endogenous mimic target to regulate the function of miRNA, and bind to miRNA through complementary sequences, thereby inhibiting the interaction between miRNA and mRNA. This inhibition of miRNA activity is called mimic target. In Arabidopsis thaliana, an endogenous lncRNA (IPS1) was found and it could bind to ath-miR399^[70]. In this study, we predicted targets and endogenous target mimics of miRNA that were established in V. vinifera. These lncRNAs may certainly regulate the corresponding miRNAs, so it is better to combine lncRNAs with miRNAs research in future studies on grape resistance.

The expression pattern of grape lncRNA was related with downy mildew stress

With regards to the RNA-seq data, the read counts of a gene region depend on the gene length and sequencing depth, so the read counts of genes or transcripts should be normalized. The FPKM statistics of grape lncRNA sequencing data showed that about 2/3 of the genes were extremely low in expression, which was consistent with previous reports in other species, and the expression level of lncRNA was lower than protein-coding genes. When P. viticola invades grape leaves, it enters through stomata and then germinates^[71]. Previous studies showed successful infection on the first day of grape leaves after inoculation, and the production of resistant response to its zoospores was observed at 1 dpi after invasion of the host cell^[72]. The number of lncRNA expressed in grape leaves inoculated was consistently higher in 'Liuba-8' than in 'Pinot noir'. Besides, the leaves expressed the most lncRNAs at 1 dpi after inoculation. This indicates that the expression of grape lncRNA may be related to downy mildew resistance. Overall, during the five days of grape leaves after inoculation, it showed that the number of DEGs in 'Liuba-8' was gradually improving and more than that in 'Pinot Noir' except on the fifth day, which is consistent with the pathogen development in the leaves of both species. Thus, it could be interesting and significant to analysis the function of genes that are differentially expressed in each stage or each species in further studies.

MSTRG.12742.1 may regulate the activity of stomata

The differentially expressed genes were screened out through RNA-seq, and lncRNAs that may have positive effects on resistance to downy mildew were selected through fluorescence quantitative PCR verification. MSTRG.12742.1 were overexpressed in PN leaves by transient transformation to detect whether it inhibited the development of P. viticola. In a previous study, it was also proved that a tomato lncRNA had a certain control effect on Phytophthora infestans by transient transformation into leaves^[13]. In the study of disease resistance genes in grape, due to the long time needed to obtain stable transgenic grape plants, preliminary verification was generally carried out in other model plants or by means of instantaneous transformation of leaves^[73]. The transient experiment has some limitations, which cannot guarantee the stable expression of genes in all leaves. At the same time, long-term in vitro culture will lead to the gradual decrease of gene expression and affect the expression of disease resistance. In the later research, it is best to obtain stably expressed plants, and then carry out the detection of inoculation and disease resistance. When P. viticola invades grape leaves, it enters through stomata, so the stomatal opening of grape leaves affects grape resistance to downy mildew to a certain extent. The expression level of MSTRG.12742.1 was significantly increased after transient transfection of PN leaves, and the protein-coding gene VIT_204s0008g02671.1, which had the highest positive correlation with the expression of MSTRG.12742.1, also showed the same result, indicating that there was a certain degree of mutual regulation between these two genes. To search VIT_204s0008g02671.1 on NCBI, it revealed the code for Cryptochrome DASH. It has been reported that cryptochrome regulates a variety of life activities in plants, including stomatal opening, and plays an important role in plant disease resistance^[74]. This was also consistent with the observed staining results that lncRNA-MSTRG.12742.1 showed a blocking effect in the early stage of P. viticola infection. At the same time, the ultrastructure study of P. viticola development in 'Pinot Noir' leaves at 1, 3, 5 and 7 dpi showed that the results were basically consistent with the observation of the colonization process of P. viticola after aniline blue staining in this experiment^[75].

One problem that should be mentioned is the lack of strict non-inoculated mock up which might cause a potential false positive, though there are some research which has reported without water mock. The samples harvested at multiple time points partly reduced the false positive rate. However, in future studies, it is better to have a much stricter mock up to ensure the reliability of the experiment.

In the early stage of this study, transcriptome sequencing was performed on the leaves of Chinese wild Vitis piasezkii accession 'Liuba-8' and European grape Vitis vinifera cultivar 'Pinot Noir' at six time points of 0, 12, 24, 48, 96 and 120 h after inoculation with downy mildew. Based on transcriptome sequencing results, possible lncRNAs were screened out in this study and the characteristics of grape lncRNAs were analyzed. Then, the cis- and trans- regulatory target genes of lncRNAs were predicted, and the possible functions of lncRNAs were analyzed. Then the microRNAs that lncRNAs might target were predicted and lncRNAs play a role as microRNA endogenous target simulators. A lncRNA that might play a positive role in grape resistance to downy mildew was screened and preliminarily verified by transient overexpression.

Overall, our research analyzed a series of characteristics of grape lncRNAs, which provided a new direction for the follow-up study of grape resistance to downy mildew. The research on lncRNA may reveal the new resistance mechanism of grape to downy mildew and provide a new theoretical basis for the rapid development of grape breeding.

The data that support the findings of this study are available from the corresponding author upon reasonable request. The raw sequence data had been deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (www.ncbi.nlm.nih.gov/sra, under accession number PRJNA22951, PRJNA902490, PRJNA902432, SUB12290370).