Abstract
BACKGROUND:
Long non-coding ribonucleic acids (lncRNAs) have been linked to many important biological processes, including fruit ripening. “Fengzao” is an early-ripening bud mutant of “Kyoho” that matures nearly 30 days earlier. However, the molecular networks controlling early ripening in Fengzao are still poorly understood in comparison to those in Kyoho.
OBJECTIVE:
The aim of the study was to gain a better understanding of the regulatory role of lncRNAs in the early ripening of grape berries.
METHODS:
The RNA-sequencing approach and quantitative real-time polymerase chain reaction validation were employed to identify lncRNAs and profile their expression patterns during berry development.
RESULTS:
In total, 24726 lncRNAs were identified, of which 19699 were differentially expressed (DE-lncRNAs) between developmental stages. The target genes of these lncRNAs and their regulatory relationships were predicted. The oxidoreductase activity, plant–pathogen interaction, plant hormone signal transduction, and flavanol biosynthesis pathways were enriched in the target genes of DE-lncRNAs, and six key lncRNAs (TCONS_00221683, TCONS_00684459, TCONS_00022149, TCONS_00167247, TCONS_00258125, and TCONS_00261813) were identified that may regulate the early ripening of grape berries.
CONCLUSIONS:
The results contribute to the understanding of the role of lncRNAs in early ripening of grape berries and will provide new insights for molecular breeding of grapes.
Introduction
Fruit ripening has received considerable scientific interest because of its uniqueness in plant biology and the fact that fruit is an important part of the human diet [1]. It is a complex and highly coordinated developmental process involving a series of changes at the physiological and biochemical levels [2]. Fruit ripening is regulated by thousands of genes that control progressive softening and/or lignification of pericarp layers, accumulation of sugars, acids and pigments, and the release of volatiles [2]. A deeper understanding of the fruit-ripening process is critical to crop improvement.
The common grape vine (Vitis vinifera L.) is one of the most widely cultivated fruiting plants in the world, and it has been grown for fruit and wine production for thousands of years. At present, most widely cultivated grape varieties in China are middle–late-ripening types, among which Kyoho is a typical representative cultivar. This cultivar has multiple advantages, including large berries, high production, and adaptation to environments with high temperatures and rainy weather [3]. Fengzao is a bud mutant of Kyoho that ripens in early July in Henan Province, China, which is nearly 30 days earlier than Kyoho [4]. The phenological, physiological, and molecular differences between Fengzao and Kyoho have been investigated in previous studies [5–7]. Histological and molecular analysis have shown that their genetic background is highly uniform [5]. RNA-sequencing (RNA-seq) analysis and MicroRNA (miRNA) analysis have demonstrated that the reactive oxygen species (ROS)-related genes and some key miRNAs are expressed differently in Fengzao and Kyoho [6–7]. In addition, a number of transcriptional analyses have been performed to investigate the ripening process in grapes at different developmental stages, across different tissues, and among different cultivars [8–10]. These studies have revealed a number of genes and pathways in grape berry ripening. However, the mechanisms by which molecular networks control early ripening in Fengzao are still poorly understood in comparison to those in Kyoho.
Long non-coding RNAs (lncRNAs) are those that are longer than 200 nt and lack a region for protein coding [11]. According to their distribution in the genome, lncRNAs can be divided into three types: long intergenic non-coding RNAs (lincRNAs), long non-coding natural antisense transcripts (lncNATs), and intronic RNAs [12]. They are involved in chromatin modification, epigenetic regulation, genomic imprinting, transcriptional control, and pre- and post-translational messenger RNA (mRNA) processing [13–14]. A vast number of lncRNAs have been identified and characterized in plants [15–16]. Functional studies have shown that lncRNAs play widespread roles in diverse biological processes, including fruit ripening, regulation of flowering time, and abiotic stress in plants [17–18]. It has been shown from RNA-sequencing analysis that many lncRNAs are significantly up- or down-regulated in tomato fruits of ripening mutants [16]. The silencing of two tomato intergenic lncRNAs resulted in an obvious delay in the ripening of wild-type fruit [19]. This indicates that lncRNAs might be essential regulators of tomato fruit ripening. The regulatory role of lncRNAs in fruit ripening in other plants has also been reported, such as Cucumis melo [20], Hippophae rhamnoides [21], and Fragaria pentaphylla [22]. However, to our knowledge, the way in which lncRNAs function during the ripening of grapes is still largely unknown.
To gain a better understanding of the regulatory role of lncRNAs in early-ripening grape berry, RNA-seq profiles taken during berry development in Kyoho and Fengzao were investigated. Differentially expressed lncRNAs during early-stage berry development were identified, and their corresponding target genes were predicted. The enrichment of the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of the target genes and the differentially expressed mRNAs were also analyzed. The results strongly suggested that lncRNAs play an important role in grape berry ripening. The findings provide new insight into the understanding of grape berry development.
Materials and methods
Plant materials
Both Fengzao and Kyoho vines with the same age were hedgerow cultivated in the same vineyard with the same viticulture management practices on the farm of Henan University of Science and Technology, Luoyang, China. The berries of the Fengzao and Kyoho vines were collected at several developmental stages from five individual vines (grafting on the same rootstock), and ten berries were taken from each vine. A total of 30 berries were collected mixed as one replicate and immediately frozen in liquid nitrogen, then stored at –80°C until further use. Based on comprehensive consideration of the grape growth stages [23] and our previous results [4, 24], the berries were sampled at the E-L 33 (hard green berries), E-L 34 (starting to soften), E-L 35 (v
RNA isolation, library construction, and sequencing
Total RNA was isolated from the berries according to the method of Rienth et al. [25]. The integrity of RNA was assessed with the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total of 3μg RNA was used for preparation of each sample. Ribosomal RNA was removed by the Ribo-Zero™ rRNA Removal Kit (Illumina, USA). Sequencing libraries were generated using the rRNA-depleted RNA by NEBNext®Ultra™ Directional RNA Library Prep Kit for Illumina® (New England Biolabs, MA, USA) following the manufacturer’s recommendations. To preferentially select cDNA fragments of 150 to 200 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, CA, USA). Then, USER Enzyme (New England Biolabs, MA, USA) was used with adaptor-ligated cDNA before polymerase chain reaction (PCR) analysis. The PCR products were purified using the AMPure XP system and library quality was assessed using the Agilent Bio analyzer 2100 system. The PCR was performed with Phusion High-Fidelity DNA Polymerase, universal PCR primers, and index primers. The libraries were sequenced on an Illumina HiSeq 2000 platform (Illumina, San Diego, CA, USA) and 100 bp paired ends were generated.
Quality control, mapping, and transcriptome assembly
Raw reads in FASTQ format were first processed through in-house Perl scripts to remove reads containing more than 10% of undetermined bases (denoted as “N”). Then, sequencing adapters and low-quality reads were removed with Trimmomatic version 0.35 [26]. The FastQC tool (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) was used to check the quality of trimmed reads, and all the subsequent analyses were based on these high-quality clean reads. The grape reference genome [27] and gene annotation files (v2.1) were downloaded from the grape genome database (http://www.genoscope.cns.fr). An index of the reference genome was built using Bowtie v2.0.6 [28] and all clean reads from each sample were aligned to the reference genome using HISAT2 v2.0.4 with the parameter “–rna-strandness RF” [29]. The distribution of known gene types was analyzed using HTSeq [30]. The mapped reads of each sample were assembled using Cufflinks v2.1.1 [31] using a reference-based approach with default parameters except “min-frags-per-transfrag = 0” and “–library-type”.
Identification and filtering of lncRNA candidates
The above assembled transcripts for all samples were combined using the CuffCompare software [31]. Based on the structural and functional characteristics of lncRNA, we applied five filtering steps to identify lncRNA candidates from all the assembled transcripts. First, single-exon transcripts falling within 500 bp up- or downstream of any protein-coding gene were removed. Next, we filtered out any transcripts less than 200 bp. All resulting transcripts were then filtered based upon the expression level calculated by Cufflinks, with expression thresholds of 2 FPKM (reads per kilobase of transcript per million mapped reads) for single-exon transcripts and 0.5 FPKM for multiple-exon transcripts. The obtained lncRNA candidates were then compared with known grape lncRNAs (http://genomes.cribi.unipd.it/DATA/V2/V2.1/lncRNA/) using CuffCompare, and transcripts matched with known lncRNAs were retained and included in the final data set. Any transcripts sharing identity with known rRNAs, tRNAs, snRNAs, snoRNAs, pre-miRNAs, and pseudogenes were removed. Finally, when compared with known mRNAs using CuffCompare, three classes of lncRNAs, including lincRNAs, intronic lncRNAs, and antisense lncRNAs, were identified from the remaining transcripts according to the results of class_code (http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/index.html).
To filter the potential protein-coding transcripts, the lncRNA transcripts found above were further processed using two tools for coding-potential analysis, Coding Potential Calculator (CPC) [32] and PfamScan [33]. The CPC assesses the extent and quality of the open reading frame (ORF) in a transcript and searches the sequences using a known protein sequence database to clarify the coding and non-coding transcripts. For CPC (v0.9-r2) analysis, the National Center for Biotechnology Information eukaryotic non-redundant protein database was used, and thee-value was set to 1×10–10. PfamScan (v1.3) with default parameters was used to identify the occurrence of any of the known protein family domains deposited in the Pfam database (release 31.0, both Pfam A and Pfam B). Transcripts with a CPC score >0 or with a Pfam hit were considered to have coding potential. Transcripts predicted to have coding potential by either of the tools were filtered out. The candidate lncRNAs without coding potential identified by both tools were considered as the final set of lncRNAs. The final data set was used for the following analyses.
Quantification and differential expression analysis
The relative abundance of both candidate lncRNAs and coding genes in each sample was computed by calculating the FPKM using Cufflinks (v2.1.1) [31]. The FPKM is calculated based on the length of the fragments and read counts mapped to each fragment. Gene FPKMs were computed by summing the FPKMs of the transcripts in each gene group.
Differentially expressed lncRNAs and coding genes in comparison groups were identified using the Cuffdiff program [31]. Cuffdiff determines differential expression in digital transcripts or gene expression data using a model based on the negative binomial distribution. For biological replicates, transcripts or genes with a adjusted p value <0.05 were assigned as differentially expressed. For non-biological replicates, adjusted p value <0.05 and the absolute value of log2 (fold change) >1 were set as the thresholds for significantly differential expression.
Target gene prediction and enrichment analysis
To explore the function of lncRNAs, we predicted their cis and trans target genes. This was carried out using an algorithm that searches for potential cis-target genes that are substantially near to lncRNAs (within 10kb upstream or downstream) using a genome browser and genome annotation. The function of the identified cis-target genes was then analyzed. The trans target genes of the lncRNAs can be predicted by correlation analysis or co-expression analysis between the lncRNAs and protein-coding genes. We calculated Pearson correlation coefficients between the lncRNAs and the coding genes. Coding genes with absolute correlation coefficients >0.95 were considered as trans target genes, and the function of these genes was analyzed through functional enrichment analysis to predict the main functions of the lncRNAs.
Gene Ontology enrichment analysis of differentially expressed genes, lncRNAs, and lncRNA target genes was performed using the R package “cluster Profiler” [34]. In this analysis, GO terms with a corrected p value <0.05 were considered significantly enriched by differentially expressed genes, lncRNAs, or lncRNA target genes. Similarly, pathways with p values <0.05 were considered significantly enriched by differentially expressed genes, lncRNAs, or lncRNA target genes.
Systematic cluster analysis and co-expression network analysis
Cluster analysis was performed as described previously [6]. The WGCNA package [35–36] was used to test whether the expression of lncRNAs and mRNAs correlated with samples.
Quantitative real-time PCR
To validate the expression of the differentially expressed lncRNAs and their roles in grape berry development, ten differently expressed lncRNAs and four mRNAs were chosen for quantitative real-time PCR (qPCR) analysis. Isolation of RNA, reverse transcription, and qPCR were performed as described previously [37]. The primer sequences are listed in Table 1.
The primers of lncRNAs and mRNAs used for qPCR verification
The primers of lncRNAs and mRNAs used for qPCR verification
Statistics and sequencing quality assessment
Fruit samples were harvested at four (FZ1, FZ2, FZ3, and FZ4) and five (KY1, KY2, KY3, KY4, and KY5) developmental stages from the Fengzao and Kyoho vines, respectively. A total of 2094312898 clean reads were obtained from the sequencing libraries with an average of around 116 million reads for each library (Table 2). The cleaned reads had a Q20 percentage (percentage of sequences with sequencing error rate lower than 1%) of over 97%. The GC content of each library was between 48.49% and 53.21%. The reads of high quality were mapped to the grape reference genome. Approximately 59 million reads were uniquely aligned to the reference, accounting for 50% of the total reads. Around 76.6% of the uniquely mapped reads were assigned to known exons.
Details of the sequencing results for each library
Details of the sequencing results for each library
After processing the assembled transcripts through five filtering steps, the lncRNAs were identified based on their structure and non-coding features. Of the 24726 lncRNAs identified, 9127 were intronic lncRNAs, 13731 were lincRNAs, and 1868 were lncNATs. These three types accounted for 36.9%, 55.5%, and 7.6% of total lncRNAs, respectively. There were also 5027 transcripts of uncertain coding potential (TUCPs).
We characterized the basic genomic features of the lncRNAs and mRNAs, including the length of the transcript and the ORF, the exon number, and the expression level. The full and ORF lengths of lncRNAs were shorter than them RNAs. The majority of lncRNAs (83.86%) were 200–300 bp, 13.96% were 301–1000 bp, and only 2.19% were >1000 bp in length. In contrast, 9.01% of mRNAs were 0–300 bp, 35.06% were 301–1000 bp, and 55.94% were >1000 bp in length (Fig. 1a). The ORF length of most lncRNAs (96.21%) was shorter than 100 bp, while only 18.01% of mRNA ORFs were shorter than 100 bp (Fig. 1b). In addition, we observed significant differences in the distribution of exon numbers between lncRNAs and mRNAs, with 96.54% of lncRNAs only containing one or two exons, while mRNAs had more exons and exon numbers distributed across a wider range (Fig. 1c). We then estimated the expression level of lncRNAs and mRNAs using FPKM the results showed that the overall expression levels of lncRNAs were lower than those of mRNAs (Fig. 1d).
Characterization of genomic features of the lncRNAs and mRNAs. a, Length distribution of the transcripts; b, Open reading frame length of the transcripts; c, Exon number of the transcripts; d, FPKM distribution of the transcripts.
To better understand the overall expression patterns of the transcripts in the nine samples, clustering analysis of the expression patterns in both Fengzao and Kyoho across various developmental stages was performed using the R package “TCseq”. The TCseq package provides a unified suite for clustering analysis and visualization of time-course sequencing data. The clustering analysis revealed 20 clusters for both Fengzao and Kyoho (Fig. 2). The transcripts and corresponding clusters are shown in Table S1, and the number of transcripts in each cluster is shown in Table 3. As shown in Fig. 2, most clusters showed similar expression patterns between Fengzao and Kyoho. There were several interesting clusters, such as clusters 4, 5, 6, 8, 10, 18, and 19, which showed differential expression patterns between Fengzao and Kyoho. The expression levels of 3373 transcripts from cluster 10 (including 2075 lncRNAs, 1036 mRNAs, and 262 TUCPs), which accounted for 6.32% of the total transcripts, were higher in the FZ1 stage than in the other stages, and the expression levels of 3057 transcripts from cluster 5 (including 1515 lncRNAs, 1238 mRNAs, and 304 TUCPs), which accounted for 6.28% of the total transcripts, were higher in the FZ2 stage than in the other stages. The expression levels of transcripts from cluster 19 (2829 transcripts, including 1356 lncRNAs, 1114 mRNAs, and 359 TUCPs) and cluster 4 (9788 transcripts, including 4114 lncRNAs, 4488 mRNAs, and 1186 TUCPs) were higher in the KY2 and KY3 stages, respectively. The expression patterns of genes were therefore dramatically different in clusters 4, 5, 10, and 19 during the early developmental stages, and in clusters 6, 8, and 18 during later developmental stages (Fig. 2). Genes in clusters 4, 5, 10, and 19 may thus contribute to early fruit development and ripening in Fengzao.
Cluster analysis of the gene expression patterns in both Fengzao and Kyoho across various developmental stages. The Y axis represents deviation of the expression levels (log-transformed) of a gene at different developmental stages from the mean expression of the gene in consideration. The X-axis represents the sampling time points.
The number of transcripts in each cluster
Both GO and KEGG enrichment analyses were performed in clusters 4, 5, 10, and 19. The GO term trihydroxystilbene synthase activity (GO:0050350) was enriched in cluster 10; 1,3-beta-D-glucan synthase activity (GO:0003843) was enriched in cluster 19; iron–sulfur cluster binding (GO:0051536), trihydroxystilbene synthase activity (GO:0050350), oxidoreductase activity (GO:0016671), carbohydrate binding (GO:0030246), and electron transfer activity (GO:0009055) were enriched in cluster 5; and 12 GO terms were enriched in cluster 4, the top five of which were microtubule motor activity (GO:0003777), double-stranded DNA binding (GO:0003690), motor activity (GO:0003774), hydroquinone:oxygen oxidoreductase activity (GO:0052716), and oxidoreductase activity (GO:0016722) (Fig. S1a).
The KEGG enrichment analysis revealed that flavonoid biosynthesis (ko00941), circadian rhythm-plant (ko04712), phenylalanine metabolism (ko00360), and the phosphatidylinositol signaling system (ko04070) were enriched in cluster 10; plant–pathogen interaction (ko04626) and non-homologous end-joining (ko03450) were enriched in cluster 19; photosynthesis (ko00195), carbon fixation in photosynthetic organisms (ko00710), circadian rhythm-plant (ko04712), carbon metabolism (ko01200), pentose phosphate pathway (ko00030), and flavonoid biosynthesis (ko00941) were enriched in cluster 5;and 11 terms were enriched in cluster4, the top five of which were phenylpropanoid biosynthesis (ko00940), DNA replication (ko03030), starch and sucrose metabolism (ko00500), fatty-acid biosynthesis (ko00061), and biotin metabolism (ko00780) (Fig. S1b).
The Cufflinks tool was used to quantify the relative abundance of transcripts by calculating FPKM, and Cuffdiff was used to conduct the differential expression analysis. A total of 54653 expressed transcripts were evaluated, including 29927 mRNAs, 19699 lncRNAs, and 5027 TUCPs. Of these transcripts, 2421 were significantly differentially expressed between Fengzao and Kyoho (adjusted p value <0.05 and the absolute value of log2 (fold change) <1). There were 1782 differentially expressed mRNAs (DE-mRNAs) among these 2421 transcripts, including 415 and 1367 mRNAs that were up- and down-regulated in Fengzao compared to Kyoho, respectively. Again, GO and KEGG enrichment analyses were performed to annotate these DE-mRNAs.
The GO analysis demonstrated that DE-mRNAs were enriched in microtubule motor activity (GO:0003777), flavonoid 3′,5′-hydroxylase (F3′5′H) activity (GO:0033772), hydrolase activity acting on glycosyl bonds (GO:0016798), hydrolase activity (GO:0016787), hydrolyzing O-glycosyl compounds (GO:0004553), motor activity (GO:0003774), oxidoreductase activity of oxygen (GO:0016709), and oxidoreductase activity (GO:0016705) (Fig. 3a). The KEGG enrichment analysis revealed that flavone and flavanol biosynthesis, galactose metabolism, DNA replication, phenylpropanoid biosynthesis, plant hormone signal transduction, flavonoid biosynthesis, starch and sucrose metabolism, and tropane, piperidine, and pyridine alkaloid biosynthesis were significantly enriched in DE-mRNAs between Fengzao and Kyoho (Fig. 3d, q value <0.1).
The GO and KEGG enrichment analysis of DE-mRNAs and DE-lncRNAs targets in comparisons between two cultivars across developmental stages. Enrichment term is represented by colored dots (blue indicates high enrichment and red indicates low enrichment). a–c, GO terms in molecular categories from GO enrichment analysis of DE-mRNAs, cis targets of DE-lncRNAs and trans targets of DE-lncRNAs, respectively; d–f, KEGG terms of DE-mRNAs, cis targets of DE-lncRNAs and trans targets of DE-lncRNAs, respectively.
To further explore the genes contributing to early ripening in Fengzao, comparisons between the two cultivars across the developmental stages before ripening were made, and these comparisons included FZ1 vs KY1, FZ1 vs KY2, FZ2 vs KY2, FZ2 vs KY3, and FZ3 vs KY3. The number of differentially expressed mRNAs was 888 in FZ1 vs KY1, 1857 in FZ1 vs KY2, 1494 in FZ2 vs KY2, 3580 in FZ2 vs KY3, and 4833 in FZ3 vs KY3. Compared with KY1, 536 mRNAs were up-regulated and 352 mRNAs were down-regulated in FZ1. Compared with KY2, 1132 mRNAs were up-regulated and 725 mRNAs were down-regulated in FZ1, and 1010 mRNAs were up-regulated and 484 mRNAs were down-regulated in FZ2. Compared with KY3, 1255 mRNAs were up-regulated and 2325 mRNAs were down-regulated in FZ2, and 1566 mRNAs were up-regulated and 3267 mRNAs were down-regulated in FZ3 (adjusted p value <0.05, Table 4). Once more, GO and KEGG enrichment analyses were performed to annotate these DE-mRNAs in each comparison.
The number of transcripts in each comparison
The GO analysis demonstrated that xyloglucan:xyloglucosyl transferase activity, pigment binding, glucosyltransferase activity, and chlorophyll binding were enriched in all three comparisons of FZ1 vs KY1, FZ1 vs KY2, and FZ2 vs KY2. The DE-mRNAs in FZ2 vs KY3 and FZ3 vs KY3 were enriched in microtubule motor activity, hydrolase activity, hydrolyzing O-glycosyl compounds, hydrolase activity acting on glycosyl bonds, and enzyme inhibitor activity, and DE-mRNAs in FZ3 vs KY3 were also enriched in pigment binding and chlorophyll binding (Fig. 3a).
The KEGG enrichment analysis revealed that protein processing in endoplasmic reticulum, photosynthesis-antenna proteins, photosynthesis, glutathione metabolism, and carbon fixation in photosynthetic organisms were enriched in at least two comparisons of FZ1 vs KY1, FZ1 vs KY2, and FZ2 vs KY2. Furthermore, seven pathways were enriched in FZ2 vs KY3 and FZ3 vs KY3, the top five of which were starch and sucrose metabolism, plant hormone signal transduction, fatty-acid degradation, galactose metabolism, and biotin metabolism (Fig. 3d).
From the significant difference analysis of the 19699 lncRNAs, there were 415 differentially expressed lncRNAs (DE-lncRNAs) in FZ vs KY, including 160 and 255 lncRNAs that were up- and down-regulated in Fengzao compared to Kyoho, respectively. In the comparisons of the early ripening stages between Fengzao and Kyoho, there were 168, 290, 208, 370, and 403 DE-lncRNAs in FZ1 vs KY1, FZ1 vs KY2, FZ2 vs KY2, FZ2 vs KY3, and FZ3 vs KY3, respectively. Compared with KY1, 101 of the 168 lncRNAs were up-regulated and the other 67 lncRNAs were down-regulated in FZ1. Compared with KY2, 128 of the 290 lncRNAs were up-regulated and the other 162 lncRNAs were down-regulated in FZ1, and 85 of the 208 lncRNAs were up-regulated and the other 123 lncRNAs were down-regulated in FZ2. Compared with KY3, 70 of the 370 lncRNAs were up-regulated and the other 300 lncRNAs were down-regulated in FZ2, and 107 of the 403 lncRNAs were up-regulated and the other 296 lncRNAs were down-regulated in FZ3 (Table 4 and Table S2).
The regulatory roles of lncRNAs on target genes were mediated by cis-and trans-acting mechanisms. To investigate the functions of lncRNAs, we first predicted the putative cis- and trans-regulatory target genes of lncRNAs. The cis targets of lncRNAs were sought within 10 kb upstream and downstream of the lncRNAs. The trans target genes of lncRNAs were predicted by correlation analysis or co-expression analysis between the lncRNAs and mRNAs. A total of 36689 cis-acting lncRNA–mRNA pairs were identified for the 10-kb flanking regions, including 16625 lncRNAs and 18062 mRNAs (Table S3). For the trans targets of lncRNAs, a total of 1131702 regulation pairs were identified, with 12579 lncRNAs and 4846 mRNAs (Table S4). Of the 20213 target mRNAs (including cis and trans) of the DE-lncRNAs, 5331 were overlapping with the DE-mRNAs (6990 in total), accounting for 76.27% of the total DE-mRNA.
To further predict the functions of DE-lncRNAs during early grape ripening, GO and KEGG analyses were performed using the cis and trans target genes of DE-lncRNAs between the two cultivars. Analysis of cis lncRNA targets revealed that there were no GO or KEGG terms significantly enriched in the comparison between Fengzao and Kyoho (corrected p value <0.05). Furthermore, we clustered the early developmental stages for Fengzao and Kyoho based on the GO (Fig. 3b) and KEGG (Fig. 3e) terms for cis-target genes of DE-lncRNAs. It is noteworthy that the GO terms glycerol channel activity, glycerol transmembrane transporter activity, and water channel activity were enriched in FZ2 vs KY2, and the GO term symporter activity was enriched in FZ3 vs KY3. From Fig. 3E, the KEGG terms glyoxylate and dicarboxylate metabolism, and amino sugar and nucleotide sugar metabolism, were enriched in FZ1 vs KY1 and FZ1 vs KY2, respectively.
For trans lncRNA targets, in the comparison between Fengzao and Kyoho, functional analysis illustrated that 12 GO terms were enriched, the top five of which were oxidoreductase activity, monooxygenase activity, heme binding, hydroquinone:oxygen oxidoreductase activity, and microtubule motor activity (Fig. 3c). We also performed the same analysis for the comparison between the two cultivars across the developmental stages before ripening. From this, there were two terms (motor activity and microtubule motor activity) that were enriched in FZ2 vs KY2, FZ2 vs KY3, and FZ3 vs KY3, and the other five terms (hydroquinone:oxygen oxidoreductase activity, monooxygenase activity, nutrient reservoir activity, oxidoreductase activity, and oxidoreductase activity, oxidizing metal ions) were enriched in FZ2 vs KY3 and FZ3 vs KY3 (Fig. 3c).
The KEGG analysis of trans lncRNA targets revealed that in the comparison between Fengzao and Kyoho there were 12 significant enriched KEGG terms, the top five of which were phenylpropanoid biosynthesis, starch and sucrose metabolism, flavone and flavonol biosynthesis, cyanoamino acid metabolism, and DNA replication. Furthermore, the DNA replication, phenylpropanoid biosynthesis, and starch and sucrose metabolism KEGG terms were enriched in the early-ripening-stage comparisons (Fig. 3f). These results indicate that pathways including starch and sucrose metabolism, flavone and flavonol biosynthesis, DNA replication, and phenylpropanoid biosynthesis are involved in the ripening of grapes.
Co-expression network analysis of genes in grape berry development
To obtain a global view of the regulatory mechanisms and the key genes related to the early ripening of grape berries, mRNAs and lncRNAs were used to construct a gene co-expression network with the WGCNA package. In total, 29 co-expression modules were developed, and the correlation coefficient of each co-expression module with the samples was calculated (Fig. 4). The modules with the highest correlation coefficients with FZ1, FZ2, FZ3, KY1, KY2, and KY3 were MElightgreen (r = 0.85), MEpaleturquoise and MEsteelblue (r = 0.73), MEcyan (r = 0.67), MEviolet (r = 0.73), MEgrey60 (r = 0.88), and MEturquoise (r = 0.98), respectively.
The relationships analysis between modules and samples by WGCNA (weighted gene co-expression network analysis). The numbers are module-sample weight correlation and corresponding p-values (in parentheses). The left panel shows the 29 modules. The color scale on the right shows module-sample correlation from –1 (green) to 1 (red).
To gain a further understanding of the seven positive correlation modules, GO and KEGG analyses of the mRNAs and target genes of lncRNAs in these modules were conducted. The results showed that terms were enriched mainly in MEsteelblue, and the terms electron transfer activity and chlorophyll binding were enriched in GO analysis (Fig. 5a), while carbon fixation in photosynthetic organisms, carbon metabolism, photosynthesis, photosynthesis-antenna proteins, fructose and mannose metabolism, and the pentose phosphate pathway were enriched in KEGG analysis (Fig. 5b). These indicate that the FZ2 stage may play an important role in the early maturing of Fengzao. To identify the genes that play pivotal roles in essential biological processes, genes in the seven modules above were analyzed by connected sub networks (Fig. 6). We defined those with more than five edges as “hub genes”, and 31 hub genes (nine mRNAs and 22 lncRNAs) were screened out, including WRKY51 (VIT_07s0031g01710), alternative oxidase gene AOX1D (VIT_02s0033g01380), and glyceraldehyde-3-phosphate dehydrogenase gene GAPA-2 (VIT_14s0068g00680). Among these, GAPA-2 is a key enzyme gene in the glycolysis pathway, which is thought to play a role in the early ripening of pears [38].
The GO and KEGG enrichment analysis of genes related to early-ripening of grape berry from co-expression network analysis. Enrichment term is represented by colored dots (blue indicates high enrichment and red indicates low enrichment). a, GO terms in molecular categories from GO enrichment analysis; b, KEGG terms.
Subnetwork analyses of important module genes in early development stages of grape berry. a, FZ1-lightgreen; b, FZ2-paleturquoise; c, FZ2-steelblue; d, FZ3-cyan; e, KY1-violet; f, KY2-grey60; g, KY3-turquoise.
The lncRNAs that appeared in all three analyses, the DE-lncRNAs (FZ vs KY, FZ1 vs KY1, FZ1 vs KY2, FZ2 vs KY2, FZ2 vs KY3, and FZ3 vs KY3), TCseq analysis (clusters 4, 5, 10, and 19), and WGCNA analysis (modules MElightgreen, MEpaleturquoise, MEsteelblue, MEcyan, MEviolet, MEgrey60, and ME turquoise), were considered to be important lncRNAs in the early developmental stages of the grape berry. The Venn diagram in Fig. 7a shows that there are six overlapping lncRNAs: TCONS_00221683, TCONS_00684459, TCONS_00022149, TCONS_00167247, TCONS_00258125, and TCONS_00261813. These six lncRNAs are considered to play key roles in the early ripening of grape berries. Similar analysis was performed to obtain key mRNAs, and there were 86 overlapping mRNAs (Table S5) in the target genes of DE-lncRNAs, DE-mRNAs, TCseq analysis, and WGCNA analysis (Fig. 7b). Caffeate O-methyltransferase activity, monooxygenase activity, aspartic-type endopeptidase activity, and aspartic-type peptidase activity were significantly enriched in the GO analysis of the 86 mRNAs (Fig. 8), but none of the KEGG terms was enriched. By comparing the 86 overlapping mRNAs and hub genes in WGCNA, it was found that three mRNAs were present in both, and all were target genes of TCONS_00684459. These are WRKY1 (VIT_07s0031g01710), alternative oxidase gene AOX1D (VIT_02s0033g01380), and calcium-binding protein KIC-like gene CBP(VIT_15s0046g00070). These suggest that TCONS_00684459 may regulate the early ripening of grape berries by regulating the important target genes.
Venn diagram showing the lncRNAs and mRNAs in this study among different analyses. a, lncRNAs among differential expression analysis, TCseq cluster analysis and WGCNA analysis; b, mRNAs between DE-mRNAs and targets of lncRNAs from DE-lncRNA, TCseq cluster analysis and WGCNA analysis.
GO terms in molecular categories from GO enrichment analysis of 86 overlapping mRNAs. Enrichment term is represented by colored dots (blue indicates high enrichment and red indicates low enrichment).
To assess the repeatability of the sequencing data, ten lncRNAs and four mRNAs were selected to determine their expression levels using qPCR. Six of the ten lncRNAs were those that were overlapping among DE-lncRNAs, TCseq analysis, and WGCNA analysis, and the remaining four were randomly selected from the DE-lncRNAs. The four mRNAs were GAPA-2 (VIT_14s0068g00680), WRKY51 (VIT_07s0031g01710), AOX1D (VIT_02s0033g01380), and CBP (VIT_15s0046g00070).
The validation tests were carried out on berry samples that were collected at the same developmental stages as those used in the RNA-seq experiment. Most tested lncRNAs showed similar expression trends in the qPCR data to those revealed by lncRNA-seq (Fig. 9). The expression of TCONS_00221683, TCONS_00022149, TCONS_00167247, and TCONS_00258125 in KY2 was higher than that in other stages in Kyoho and all stages detected in Fengzao, and the expression of TCONS_00261813 was the highest in KY3. The expression of TCONS_00684459 and its target genes (VIT_07s0031g01710, VIT_02s0033g01380, and VIT_15s0046g00070) was higher in FZ1 than in KY1. This result illustrates that our high-throughput data were reliable.
qPCR validation of lncRNAs and mRNAs at different berry development stages in Fengzao and Kyoho. Relative quantity is based on the expression of the reference gene ubiquitin 1. X-axis indicates different stages and Y-axis indicates the expression of lncRNAs or mRNAs relative to reference gene. Data are mean±SD from three biological replicates. *, P < 0.05; **, P < 0.01 by student t test.
Recently, lncRNAs have been recognized as an important class of regulatory molecules. Plant lncRNAs are linked to biological processes such as gene silencing, regulation of flowering time, abiotic stress response, and fruit developmental pathways [39]. In previous studies, lncRNAs related to fruit development have been identified in tomato [16, 19], muskmelon [20], and strawberry [22], among others. However, the lncRNAs involved in the grape berry ripening process have remained largely unknown. To characterize the lncRNAs relevant to grape berry development, lncRNAs and their expression profiles across the key developmental stages of Kyoho and Fengzao were examined. A total of 24726 lncRNAs were identified, which included 9127 intronic lncRNAs, 13731 lincRNAs, and 1868 lncNATs. Furthermore, structural analysis showed that these lncRNAs were shorter than the mRNAs, and the lncRNAs expression levels were lower than those of the mRNAs (Fig. 1). The general features were consistent with previous studies [16, 40]. The TCseq clustering analysis classified lncRNAs into 20 clusters (Fig. 2), and the expression levels of lncRNAs varied at different developmental stages, indicating that lncRNAs play diverse roles during grape berry development. Interestingly, the expression levels of lncRNAs in clusters 4, 5, 10, and 19 changed dramatically during the early developmental stages, so they may contribute to early fruit development and ripening in Fengzao.
In our previous study, RNA-seq and miRNAs analyses were carried out on Fengzao and Kyoho during berry development, and the results indicated that superoxide dismutase (SOD) plays an important role in the early ripening of grape berries. The activity of the SOD enzyme in Fengzao was lower than that in Kyoho, except at the v
Pathogenesis-related (PR) proteins, which are generally considered as plant defense proteins, are believed to be involved in fruit ripening. Several differentially expressed genes that are related to disease-resistance proteins were observed in Fengzao and Kyoho during berry development, and the disease-resistance protein genes VIT_213s0067g01100 and tobacco mosaic virus resistance protein N-like gene VIT_200s0238g00060 were among the most differentially expressed genes detected in the study [6]. In this study, the KEGG term plant–pathogen interaction (ko04626) was enriched in cluster 19 of the TCseq analysis (Fig. S1b). Among the six overlapping lncRNAs, which appeared in the DE-lncRNAs, the TCseq analysis, and the WGCNA analysis (Fig. 7a), three (TCONS_00221683, TCONS_00022149, and TCONS_00167247) share the same target gene, VIT_03s0088g00910 (the basic form of PR protein 1), which is considered as a key gene in the early ripening of the Fengzao berry. Lijavetzky et al. [44] also observed the activation of pathogen-defense gene expression responses in the pericarp upon berry ripeness in “Muscat Hamburg”. Andthe role of PR proteins in the ripening of other fruits has also been demonstrated [45–46].
The only miRNA to be differentially expressed during all berry-development stages of Fengzao and Kyoho was miR159a [7]. The grape targets of miR159a include the GAMYB transcription factors MYB33, MYB65, and MYB101 [47], which participate in the signaling process induced by abscisic acid (ABA) accumulation in the presence of stress [48]. This suggests that the effect of miR159a on the development of grape berries is realized by a regulating-hormone signal pathway. Plant hormones are considered to be closely linked with fruit development and ripening [49]. Through the response of different hormones, the significant ripening regulations seem to be controlled primarily by ethylene, ABA, and auxin [50–52]. In this study, the KEGG term plant hormone signal transduction was enriched in DE-mRNAs, and the trans target genes of DE-lncRNAs were compared between two cultivars (Fig. 3d and 3f). There were 28, 57, and 80 DE-mRNAs that belonged to plant hormone signal transduction in the comparisons FZ vs KY, FZ2 vs KY3, and FZ3 vs KY3, respectively (Fig. 3d). Among the trans target genes of the DE-lncRNAs in FZ vs KY and FZ2 vs KY3, there were 41 and 25 target genes that belonged to plant hormone signal transduction, respectively (Fig. 3F). A similar result was observed in Cucumis melo, in which the target genes of DE-lncRNAs involved in fruit ripening were enriched plant hormone signal transduction terms [20], which is consistent with our results.
Other pathways have been demonstrated to be involved in fruit-ripening processes, for example cell-wall modification, oxidative stress, photosynthesis and glycolysis in pears [38], metabolic pathways, spliceosome, and biosynthesis of phenylpropanoids in oranges [53]. In this study, the DE-mRNAs and the target genes of DE-lncRNAs between Kyoho and its early-ripening bud Fengzao were found to be enriched in F3′5′H activity, flavanol biosynthesis, phenylpropanoid biosynthesis, and starch and sucrose metabolism pathways (Fig. 3d and Fig. 3f), which indicates that these pathways are involved in the early ripening of grape berries. The initiation of fruit ripening is mainly characterized by color change, hormone synthesis, and cell-wall dynamics [52]. Anthocyanins are secondary metabolites that play a substantial role in pigmentation and have many health benefits as antioxidants, including anti-tumor properties [54], and they accumulate in the early stage of fruit ripening [55–56]. In this study, three pathways related to anthocyanin synthesis (F3′5′H activity, flavanol biosynthesis, and phenylpropanoid biosynthesis) were enriched in DE-mRNAs and the target genes of DE-lncRNAs between Kyoho and Fengzao (Fig. 3d and 3f). These results showed that anthocyanin synthesis plays an important role in the early ripening of Fengzao, which is consistent with previous research about apple fruit development. Three different stages of apple development were used for transcriptome assembly using RNA-seq, and differentially expressed genes were involved in hormonal signaling pathways and anthocyanin biosynthesis pathways [57].
Combining differential expression analysis, TCseq analysis, and WGCNA analysis, it was found that six lncRNAs and 86 mRNAs may play key roles in the early ripening of Fengzao (Fig. 7). The six lncRNAs are TCONS_00221683, TCONS_00684459, TCONS_00022149, TCONS_00167247, TCONS_00258125, and TCONS_00261813. The 86 mRNAs not only contained genes in key hormone signal transduction pathways (VIT_03s0038g01150, VIT_05s0020g01070, VIT_04s0023g02350, and VIT_03s0038g01180), anthocyanin synthesis (VIT_06s0009g03140), starch and sucrose metabolism (VIT_12s0055g00310, VIT_06s0004g04330, and VIT_07s0289g00080) mentioned earlier, but also included methyltransferase (VIT_14s0066g01450, VIT_15s0048g02490, and VIT_15s0048g02450) and WRKY transcription factors (VIT_08s0058g01390, VIT_07s0031g01710, and VIT_04s0008g01470) (Table S2).
Methylation of DNA is an epigenetic modification that affects many biological processes and is widely presented in animals and plants. Some studies have revealed that DNA methylation plays a very important role in the growth and development of fruit [58–59]. In our previous study, different concentrations of azacitidine (5-azaC), an inhibitor of methyltransferase, were sprayed onto Kyoho grape berries, and the results showed that treatment with 100μM 5-azaC could cause berry ripening in Kyohoto occur 20 days earlier than in the control [60]. In this study, three genes (VIT_14s0066g01450, VIT_15s0048g02490, and VIT_15s0048g02450) related to methylation appeared to play important roles in the early ripening of Fengzao, which further confirms the role of methylation in the development and early ripening of grape berries.
The WRKY-family genes are generally considered to be stress tolerance regulators, but recent studies also implicate WRKY transcription factors in fruit development. It has been predicted that FvWRKY genes may operate in the ABA, indole-3-acetic acid, and sucrose signaling networks during strawberry fruit development [61]. The differentially expressed WRKY70 gene was found to be co-expressed with several DE-lncRNAs in the transcriptome of Cucumis melo fruit at four developmental stages, so it was speculated that WRKY70 might be a negative regulatory factor in the process of fruit development [20]. Of the 86 important mRNAs, three are WRKY transcription factors (Table S2), which indicates that WRKY plays a role in grape berry development. Further analysis of the 86 overlapping mRNAs and hub genes in WGCNA found that WRKY1 (VIT_07s0031g01710), alternative oxidase gene AOX1D (VIT_02s0033g01380), and calcium-binding protein KIC-like gene CBP (VIT_15s0046g00070) was present in both of them, and this shows that they are key genes in the early-ripening process of Fengzao. Interestingly, WRKY1, AOX1D, and CBP are all target genes of TCONS_00684459. Results from RNA-seq and qPCR showed that the expression of TCONS_00684459, WRKY1, AOX1D, and CBP was much higher in FZ1 than in other stages (Fig. 9). These results suggest that TCONS_00684459 may regulate the early ripening of grape berries by regulating the important target genes. However, the function of TCONS_00684459 and its regulatory role on targets needs to be explored in more depth.
Conclusions
Fruit ripening is a well-regulated process in which many pathways and genes are probably involved. In this study, the oxidoreductase activity, plant–pathogen interaction, plant hormone signal transduction, flavanol biosynthesis, phenylpropanoid biosynthesis, and starch and sucrose metabolism pathways were found to be enriched in the target genes of DE-lncRNAs. Six key lncRNAs were identified (TCONS_00221683, TCONS_00684459, TCONS_00022149, TCONS_00167247, TCONS_00258125, and TCONS_00261813) that may regulate the early ripening of grape berries by regulating the important target genes. The results contribute to the understanding of the role of lncRNAs in the early ripening of grape berries and will provide new insights for molecular breeding of grapes.
Funding
This work was financially supported by Natural Science Foundation of China (NSFC: 31672106, U1904113), Zhongyuan Science and Technology Innovation Leaders (194200510007), National Key Research and Development Program of China (2018YFD1000105), Key Scientific Research Projects of Colleges and Universities in Henan Province (20A210009), and Research Funding for Young Academic Leaders of Henan University of Science and Technology (4026-13490004).
Conflict of interest
The authors have no conflict of interest to report.
Footnotes
Acknowledgments
This work was financially supported by Natural Science Foundation of China (NSFC: 31672106, U1904113), Zhongyuan Science and Technology Innovation Leaders (194200510007), National Key Research and Development Program of China (2018YFD1000105), Key Scientific Research Projects of Colleges and Universities in Henan Province (20A210009), and Research Funding for Young Academic Leaders of Henan University of Science and Technology (4026-13490004). The Funding bodies were not involved in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.