Genome-wide identification and characterization of long non-coding RNAs involved in fruit ripening and the climacteric in Cucumis melo
BMC Plant Biology volume 19, Article number: 369 (2019)
Cucumis melo is a suitable study material for investigation of fruit ripening owing to its climacteric nature. Long non-coding RNAs have been linked to many important biological processes, such as fruit ripening, flowering time regulation, and abiotic stress responses in plants. However, knowledge of the regulatory roles of lncRNAs underlying the ripening process in C. melo are largely unknown. In this study the complete transcriptome of Cucumis melo L. cv. Hetao fruit at four developmental stages was sequenced and analyzed. The potential role of lncRNAs was predicted based on the function of differentially expressed target genes and correlated genes.
In total, 3857 lncRNAs were assembled and annotated, of which 1601 were differentially expressed between developmental stages. The target genes of these lncRNAs and the regulatory relationship (cis- or trans-acting) were predicted. The target genes were enriched with GO terms for biological process, such as response to auxin stimulus and hormone biosynthetic process. Enriched KEGG pathways included plant hormone signal transduction and carotenoid biosynthesis. Co-expression network construction showed that LNC_002345 and LNC_000154, which were highly expressed, might co-regulate with mutiple genes associated with auxin signal transduction and acted in the same pathways. We identified lncRNAs (LNC_000987, LNC_000693, LNC_001323, LNC_003610, LNC_001263 and LNC_003380) that were correlated with fruit ripening and the climacteric, and may participate in the regulation of ethylene biosynthesis and metabolism and the ABA signaling pathway. A number of crucial transcription factors, such as ERFs, WRKY70, NAC56, and NAC72, may also play important roles in the regulation of fruit ripening in C. melo.
Our results predict the regulatory functions of the lncRNAs during melon fruit development and ripening, and 142 highly expressed lncRNAs (average FPKM > 100) were identified. These lncRNAs participate in the regulation of auxin signal transduction, ethylene, sucrose biosynthesis and metabolism, the ABA signaling pathway, and transcription factors, thus regulating fruit development and ripening.
Genome-wide transcriptome sequencing has revealed that almost 90% of eukaryotic genomes can be transcribed , nevertheless only 1–2% of the genome encodes proteins . Thus, non-coding RNAs (ncRNAs) constitute a dominant proportion of the transcriptome. In contrast to protein-coding mRNAs, ncRNAs are characterized by a low level of expression and sequence conservation, thus ncRNAs were originally considered as transcriptional “noise” . With the development of high-throughput sequencing technology and bioinformatic approaches, an increasing number of ncRNAs have been identified and characterized. Two groups of ncRNAs are distinguished on the basis of their length, namely small ncRNAs (sncRNAs) and long ncRNAs (lncRNAs). The sncRNAs are fewer than 200 nucleotides in length and consist of microRNAs and small nucleolar RNAs, and have been well studied in the last decade. In contrast, lncRNAs are longer than 200 nucleotides but lack protein-coding potential .
In recent years, it has been reported that lncRNAs may be transcribed from any position of a genome by RNA polymerase II, IV, and V [4,5,6]. Based on the location and context in the genome, lncRNAs can be classified into five groups, comprising sense lncRNAs, introns originating from natural antisense transcripts, long intergenic non-coding RNAs, intronic ncRNAs, and bidirectional lncRNAs . Previous reports have shown that lncRNAs play important roles in several biological processes, such as gene transcription , post-transcriptional modification , translation , transcriptional interference, protein modification, and regulation of DNA methylation [10, 11]. The best-characterized lncRNAs are mainly from humans and animals, whereas research on plant lncRNAs lags far behind in comparison. Fortunately, the accumulation of plant genetic resources in public databases has stimulated studies of plant lncRNAs. A vast number of lncRNAs have been identified and characterized in plants, such as Arabidopsis thaliana , Oryza sativa , Zea mays , and Solanum lycopersicum . For example, about 6480 A. thaliana lncRNAs were identified from 200 organ-specific and stress-response transcriptomic data sets . Functional studies have shown that lncRNAs regulate plant growth, development, reproduction, and stress responses . Swiezewski et al. reported a series of cold-induced long antisense RNAs that play a role in the epigenetic silencing of the flowering-time regulatory gene FLOWERING LOCUS C in A. thaliana . The WRKY1-activated lncRNA33732 enhances broad-spectrum resistance to pathogens in tomato . However, knowledge of plant lncRNAs is limited with respect to the diversity and distribution in the genome. Therefore, investigations of lncRNAs on additional plant species that express representative metabolic pathways and/or characteristics are required.
Cucumis melo, a member of the Cucurbitaceae family, is an attractive model plant for investigation of important biological processes, including fruit ripening, sex determination, and stress tolerance [18,19,20]. Cucumis melo fruit show notable variation in ripening physiology. The fruit of different varieties can be categorized as climacteric or non-climacteric based on the ripening-related respiration rate . The climacteric is the final physiological stage that marks the beginning of climacteric-dependent fruit ripening, and results in changes in diverse internal and external traits, such as flesh softening, aroma, abscission, rind color [22, 23]. An abrupt increase in ethylene production in the fruit is characteristic of climacteric fruit ripening and usually occurs without an external influence . The physiology and molecular mechanisms of melon fruit ripening have been a research focus in recent years . Guo et al. reported that the rind color changed to yellow from green and flesh firmness increased initially, then declined significantly from 25 days after anthesis (DAA). Coincident with these changes, contents of soluble solids, sucrose, alcohols, and acids were increased in the flesh . Transcriptome profiling of Hami melon fruit development identified 7892 differentially expressed mRNAs (DE-mRNAs), including genes associated with hormone stimulus response, ethylene and sucrose biosynthesis (e.g. sucrose synthases [SUS] and 1-amino-cyclopropane-1-carboxylate oxidase [ACO]), and several transcription factor families . In addition, our laboratory reported that silencing CmACO1 postponed fruit ripening of C. melo . However, knowledge of the molecular mechanisms underlying melon fruit ripening remains extremely limited.
The melon cultivar used for this work was Cucumis melo L. cv. Hetao, white flesh, the peel turned yellow when ripe. With consideration of the hypothesis that lncRNAs might play important roles in melon fruit ripening, in this study we performed RNA sequencing (RNA-seq) to investigate lncRNAs involved in melon fruit development, and predicted the function of lncRNAs based on the position or expression patterns with their target genes. Differential gene expression was characterized among four developmental stages to gain additional insight into the genetic regulation of fruit development. Our results predict the roles of candidate genes and regulatory lncRNAs in fruit development and ripening, and provide a foundation for future studies of the molecular mechanism of ripening in climacteric fruit of melon.
Phenotypic data of the C. melo fruits
The fruit pressure test result showed that the fruits turn soft from stage G to P (Fig. 1a). Soluble solids content progressively increased from stage G to C (Fig. 1b). And the respiration rate curve peaked at 0.25–0.28 CO2 ppm/min/g (20 h and 22 h after harvest, Fig. 1c).
RNA sequencing output and genome mapping
In total, 96,820,869 paired-end reads with a mean length of 125 bp were generated in 12 libraries. After filtering out adaptor sequences and low-quality reads, 92,957,247 clean reads were obtained per library. The Q20 and Q30 of the clean reads were greater than 90%, and GC contents of the sequencing outputs were 43–45%. The percentage of clean reads in each library ranged from 92.69 to 96.87%, and approximately 66.39–77.88% of the reads were mapped to the melon reference genome (Table 1).
Characterization of the candidate lncRNAs
By matching all of the transcripts using the Coding Potential Calculator (CPC) and the Pfam database, a total of 3857 candidate lncRNAs were predicted (Fig. 2a), which comprised 3307 (85.7%) intergenic lncRNAs and 550 (14.3%) anti-sense lncRNAs (Nature anti-sense transcripts, NATs) (Fig. 2b, Additional file 1). Among the predicted lncRNAs, 142 with mean fragments per kilobase of exons per million mapped fragments (FPKM) > 100 were highly expressed. The overall length of the lncRNAs and their open reading frames were shorter than those of the mRNAs, and the exon number of the lncRNAs was less than that of the mRNAs (Fig. 2d).
Identification of differentially expressed lncRNAs
On the basis of FPKM values, a total of 1601 differentially expressed lncRNAs (DE-lncRNAs) were identified. The number of DE-lncRNAs at each developmental stage was 528 (G-vs-R), 510 (R-vs-C), 76 (C-vs-P), and 997 (G-vs-C) (Fig. 2c). The cluster analysis showed that the replicates for each developmental stage clustered together. Replicates for the R, C, and P stages formed one group, while G stage replicates formed a separate group (Additional file 11: Figure S1).
Enrichment analysis of cis-acting lncRNAs
After mapping the lncRNAs to the melon reference genome, 3854 lncRNAs were identified close to (< 100 kb) 18,277 protein-coding genes (Additional file 2), named as cis-acting lncRNAs. Gene Ontology (GO) enrichment analysis showed that seven, eight, and four GO terms were significantly enriched between stages G and C, G and R, and G and P.
All of the enriched GO terms between the G and C stages were classified under biological process, including response to stimulus (613 genes), response to hormone stimulus (50 genes), response to endogenous stimulus (50 genes), response to auxin stimulus (45 genes), hormone biosynthetic process (15 genes), pheromone metabolic process (15 genes), and pheromone biosynthetic process (15 genes; Additional file 3). The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that the target genes of lncRNAs were enriched in 119 KEGG pathways. Enriched pathways associated with fruit development included plant hormone signal transduction, carotenoid biosynthesis, carbon metabolism, and carbon fixation in photosynthetic organisms (Additional file 4).
Enriched cis-target genes in plant hormone signal transduction
A total of 171 DE-mRNAs were involved in plant hormone signal transduction. Of these DE-mRNAs, 139 co-located with 250 lncRNAs (Additional file 5), and 109 of the lncRNAs were differentially expressed. Among these DE-lncRNAs, LNC_003452 was the most highly expressed and co-located with 43 mRNAs. CmACCO3 was a target mRNA of LNC_003452, and located 94,591 bp upstream of the lncRNA.
Enrichment analysis of co-expressed genes of lncRNAs
The potential correlated pairs of lncRNAs and mRNAs were predicted using co-expression analysis. A total of 245,368 interactions between 2258 lncRNAs and 11,102 protein-coding transcripts were detected. Among all the interactions, 245,100 were trans-acting (distance beyond 100 kb, Additional file 2). Functional analysis showed that the co-expressed genes were enriched in 3788 GO terms. Among the GO terms, 2127, 1170, and 491 were classified under biological process, molecular function, and cellular component, respectively (Additional file 6). Interestingly, some of the GO terms were closely associated with fruit development, including response to hormone stimulus, response to auxin stimulus, and hormone biosynthetic process. The co-expressed genes were enriched in 120 KEGG pathways (Top 20 enriched pathways of each group showed in Additional file 7). Among the DE-mRNAs in any pairwise comparison, 171 were involved in plant hormone signal transduction pathways, of which 130 were co-expressed with 422 DE-lncRNAs. In the carotenoid biosynthesis pathway, 18 DE-mRNAs were co-expressed with 204 DE-lncRNAs.
Systematic cluster analysis, co-expression network construction, and PPI analysis
To gain insight into the temporal and spatial transcription patterns and putative functions between the lncRNAs and their co-expressed mRNAs during C. melo fruit development, 83 target genes, 91 DE-TFs (fold change > 2), and their co-expressed DE-lncRNAs were selected and a systematic cluster analysis was performed. Interestingly, the expression patterns were consistent and for each data set two groups were resolved (stages G and R, and stages C and P; Figs. 3 and 4), which indicates a correlation with development and the respiratory climacteric process in C. melo fruit. Network analysis showed that one lncRNA could co-regulated with multiple mRNAs and vice versa, suggesting complex interactions between lncRNAs and mRNAs. The protein interactions (PPI) analysis showed that CmACO1 and Sucrose Synthase Y (SUSY) interacted with ACS1 and Agamous-like 2 (AGL2), respectively. Serine/threonine-protein kinase 3 (CmSAPK3) and CmSAPK4 both interacted with CmACC1.
Expression profiling of the key elements associated with Cucumis melo fruit ripening
Genes that play important roles in pathways including hormone stimulus, carotenoid or ethylene biosynthesis, and sugar and main organic acid metabolism were selected and their expression profiles were analyzed. CmACO1, MAOX (similar to NADP-dependent malic enzyme), EBF1, CmACO7, ETR2, and SWT3B had the highest expression levels at stage C. The average FPKM of CmACO1 was extremely high (almost 13,613), followed by MAOX, EBF1, CmACO7, ETR2, and SWT3B, with average FPKM values of 1323.5, 1082, 965.77, 457.7, and 416, respectively (Fig. 5). LNC_000987, LNC_001323 and LNC_003610showed similar expression patterns at stage C, and their average FPKM values were 1954.4, 901.0 and 563.0 respectively. LNC_000693 and LNC_003380 were also highly expressed at stage C, and increased from stage G to P (Fig. 6). IAA14, LNC_002345, LNC_000154, LNC_003726, and LNC_000126 showed the highest expression levels at stage G, and had mean FPKM values of 954, 2670.6, 2533.8, 2266, and 919, respectively. SUS2 and SUSY (sucrose synthase) were up-regulated from stage G to P, and 11 DE-lncRNAs were co-expressed with SUS2. The correlations between the DE-mRNAs and lncRNAs described above are shown in the co-expression network (Fig. 7).
Among the analyzed DE-TFs, bHLH130 was significantly up-regulated from stage R to C (fold change > 4.5), and LNC_003066 was co-expressed with this DE-TF. bHLH48 and bHLH93 were highly expressed at stage G, then significantly down-regulated at stage R. The average FPKM value of NAC72 was 933 at stage C. NAC56 was highly expressed at stages R, C, and P, but the co-expressed lncRNAs were all negatively correlated (Fig. 7). WRKY70 was highly expressed at stages G and R, then significantly down-regulated at stage C. Several co-expressed lncRNAs of WRKY70 were identified, and all the correlations were positive except for one (LNC_003467). Among them, LNC_001263 was down-regulated from the G to C stages, and the FPKM value decreased from 269 to 2. ERF12, ERF76, and ERF80 were significantly down-regulated from stages G to C (Fig. 5). The gene id of the genes included were provided in Additional file 8.
Validation of differentially expressed lncRNAs
The expression patterns of 10 lncRNAs with only one transcript were validated at the four developmental stages by qRT-PCR. The results confirmed that the expression patterns of the selected lncRNAs were consistent with the expression levels calculated from the RNA-seq data (Additional file 12: Figure S2). The results showed that the pipeline used in the present study was extremely strict in identifying putative lncRNAs, and suggested that the identified lncRNAs were genuinely expressed.
An informative dataset of lncRNAs associated with fruit ripening in Cucumis melo
Increasingly, lncRNAs are recognized as an important class of regulatory molecules, but they remain poorly studied in model plants. Plant lncRNAs are linked to biological processes such as gene silencing, flowering-time regulation, abiotic stress response, and numerous other developmental pathways . However, in C. melo fruit, no information on lncRNAs was available. In the current study, by applying RNA-seq and genome mapping approaches, 3857 lncRNAs were identified at four stages of C. melo fruit development, of which 1601 lncRNAs were differentially expressed. Of all pairwise comparisons among the four analyzed stages, stages G and P had the highest number of DE-lncRNAs (1184), suggesting that gene expression regulatory patterns changed significantly in the process of C. melo fruit development. In addition, 510 DE-lncRNAs between stages R and C were identified, of which some might be closely associated with molecular regulation of the respiratory climacteric. Based on the physical location and expression relationships between lncRNAs and their target mRNAs, further prediction of the cis or trans regulatory activity of the identified lncRNAs was performed. As a result, 3854 lncRNAs were identified as close to 18,277 protein-coding genes. A total of 245,368 interactions between 2258 lncRNAs and 11,102 co-expression genes, and 245,100 were trans-acting. Finally, we analyzed the functions of DE-lncRNAs based on their target genes and the co-expression genes, of which some were involved in diverse processes associated with fruit development, including response to hormone stimulus, carotenoid or ethylene biosynthesis, and sugar and main organic acid metabolism. The present results predicted the regulatory function of the lncRNAs during C. melo fruit development and ripening, and lay a foundation for further detailed studies.
Cis/trans roles of the lncRNAs and their target DE-genes involved in plant hormone signal transduction
Fruits are important resources for human and animal nutrition owing to their high vitamin, mineral, and sugar contents. In addition, fruits are rich sources of antioxidant compounds, such as carotenoids, anthocyanins, and flavonoids , which are of great economic value. Ripening is crucial for development of the flavor and nutritional quality of fruits . Fruit ripening is a developmental process coordinated by complex networks of interacting genes and signaling pathways. Plant hormones play important roles in the regulation of this process. In fruit trees and crops, a variety of genes and regulators encode and/or regulate enzymes in the plant hormone signaling transduction pathways of fruit components during ripening. For example, in tomato, the MADS-box genes ripening inhibitor (RIN) and agamous-like 1 (AGL1) dramatically affect fruit ripening . In tomato, pepper, banana, and strawberry, the dominant endogenous auxin indole-3-acetic acid (IAA) defers early fruit development . However, current knowledge of the roles of hormones (except for ethylene) in the development and ripening of climacteric and non-climacteric fruits is limited, therefore the genetic and molecular factors involved in the plant hormone signal transduction pathway that regulate ripening should be identified and characterized. One of the main objectives in the present study was to understand how the co-expression network among the lncRNAs and mRNAs controls development of melon fruit. As a result, we identified 109 DE-lncRNAs co-located with the DE-mRNAs involved in plant hormone signal transduction, and 422 DE-lncRNAs co-expressed with the DE-mRNAs involved in the same process.
Auxin plays a primary role during the conversion of the ovary into a growing fruit and controls many aspects of fruit development, including fruit set and growth, ripening and abscission . In auxin transduction, auxin-induced proteins and IAA are the most well-studied components and are encoded by Auxin-Induced (AX) genes and IAA. In the current study, 22 auxin transduction-related DE-genes were all co-expressed with LNC_002345, including IAA4, IAA8, IAA9, IAA14, AX22D, and AX15A, and 16 DE-genes were co-expressed with LNC_000154, including AX22D, AX15A, IAA8, and GID1C (a gibberellin receptor) (Fig. 6). These genes were all highly expressed at stage G, then significantly down-regulated at stage R. However, LNC_002345 was highly expressed at all four stages and the average FPKM value was 2670.6 at stage G, then was down-regulated by more than two-times at stage R and its expression remained stable until stage P. The FPKM value of LNC_000154 at stage G was 2533.8, then the gene was significantly down-regulated to 395 at stage R. We predicted that LNC_002345 and LNC_000154 might regulate the expression of genes associated with auxin signal transduction, and thereby be involved in fruit development. Other highly expressed lncRNAs that also interacted with auxin signal transduction genes are shown in Fig. 7.
With regard to the abscisic acid (ABA) signaling pathway, ABA promotes strawberry fruit ripening, and the ABA content peaks earlier than that of ethylene, which suggests that ABA and ethylene both promote strawberry fruit maturation . SAPK genes encode serine/threonine-protein kinase and play a role in the ABA signaling pathway. Among the present findings, the most highly expressed lncRNA at stage C was LNC_000693, which was co-expressed with CmSAPK3. LNC_003380 and LNC_000415 were up-regulated from stage G to stage P (average FPKM values increased from 8.7 to 849, and 0.5 to 1124, respectively), and all three lncRNAs were co-expressed with CmSAPK4. The expression patterns and interactions between lncRNAs and SAPK genes of the above-mentioned three lncRNAs suggest that they might regulate C. melo fruit ripening by participating in the ABA signaling pathway. PPI analysis showed that CmSAPK3 and CmSAPK4 both interacted with CmACC1 (which encodes an ethylene precursor), implying that SAPK genes may be associated with ethylene biosynthesis in C. melo fruit.
lncRNAs and their co-expressed genes associated with ethylene and sucrose biosynthesis and signal transduction
Fleshy fruits can be classified into two physiological groups, climacteric and non-climacteric, according to their respiration characteristics and rate of ethylene release during ripening. Ethylene synthesis in climacteric fruits, such as tomato, apple, banana, and C. melo, is essential for the normal fruit ripening process . In the present study, CmACO1 and CmACO7 were significantly up-regulated at the post-ripening stage P, which was consistent with previous reports . Mutation of ETR1–1 results in disruption of ethylene binding during the ethylene response in Arabidopsis . In C. melo, EBF1, ETR1, and ETR2 were up-regulated during fruit development. We predicted that ETR1 might interact with LNC_000866 and LNC_001504. The average FPKM value of LNC_000866 from stage G to C increased from 99.8 to 313, then was down-regulated to 168 at stage P. ETR2 might interact with LNC_003066 in the same pathway. These results revealed that lncRNAs may interact with ETR1/2 and play roles in ethylene signaling transduction in C. melo.
Recent studies have acknowledged that sucrose acts as a signal to modulate a wide range of processes in plants, including fruit development. Abdullah et al. detected four genes of the SUS family that were significantly up-regulated during fruit development of Chinese pear . In the present study, SUS2 and SUSY (which encode sucrose synthases) showed similar expression patterns from stages G to P, which was consistent with the physiological data (soluble solids content test). In addition, the expression levels of CmMAOX (NADP-dependent malic enzyme) and CmSWT3B (bidirectional sugar transporter) increased by 52 and 75 times, respectively, from stages G to C, and then were down-regulated at stage P. The changes in expression of these two genes suggest that they may perform regulatory roles during fruit development in C. melo, but the specific functions of the two genes have not been elucidated and require investigation in the near future.
Expression profiling and regulatory role of transcription factors and their co-expressed lncRNAs in Cucumis melo
Transcription factors (TFs) are regulatory proteins responsible for the transcriptional activation or repression of target genes by recognizing specific, short DNA sequences (6–12 bp) in their regulatory regions. Studies on the regulatory mechanisms and functions of TFs associated with fruit development have increasingly been considered as research hotspots. However, current knowledge is extremely limited and is mainly confined to several TF families, such as ERF, RIN, and MADS-box , thus an improved understanding of the roles of TFs in the process of fruit ripening is needed.
In the present study, we analyzed the DE-TFs during C. melo fruit ripening, and determined that the TFs differentially expressed among the four developmental stages were largely distributed in the bHLH, ERF, NAC, and WRKY TF families. The present results were consistent with observations reported for tomato fruit ripening , suggesting that the transcriptional strategies during ripening in C. melo may be similar with other climacteric fruit. The systematic cluster analysis showed that the expression patterns of DE-TFs (fold change > 2) and their co-expressed DE-lncRNAs are consistent (Fig. 4), which indicated that these DE-TFs may be associated with regulation of the respiratory climacteric in C. melo.
The NAC gene family is considered to be involved in the regulation of fruit ripening. NAC protein can bind to the promoter of ACS and ACO in ripening peach fruit tissues . NAC4 has been confirmed to play a role in tomato fruit ripening . Ríos et al. identified NAC56 as a fruit ripening quantitative trait locus that regulates climacteric ripening in melon . In the present investigation, NAC56 was significantly up-regulated from stages G to R, which was consistent with the expected results, thus the gene may play a role in regulation of fruit ripening in C. melo. In addition, NAC72 was remarkably up-regulated from stages R to C, then down-regulated slightly at stage P. We speculate that NAC72 might play a crucial role in the C. melo climacteric process. The bHLH gene family is a common TF family among plants. MdbHLH93 interacts directly with an ABA-responsive protein, MdBT2, and thus delays leaf senescence . In C. melo, bHLH130 was significantly up-regulated from stages R to C (fold change > 4.5), and LNC_003066 was co-expressed with bHLH130. LNC_003066 and bHLH130 might function as positive regulators during C. melo fruit ripening. bHLH48 and bHLH93 were highly expressed at stage G, then were down-regulated significantly at stage R; in addition, we detected that LNC_000126 was co-expressed with bHLH93. Whether LNC_000126 and bHLH93 are negative regulators of C. melo fruit ripening requires further study. WRKY family genes are generally considered to be stress tolerance regulators, but recent studies implicate WRKY TFs in fruit development. It has been confirmed that the WRKY domain facilitates binding of the protein to the SURE (sugar-responsive cis-element) element in the promoter regions of target genes . Zhou et al. predicted that FvWRKY genes may operate in the ABA, IAA, and sucrose signaling network during strawberry fruit development . We determined that the average FPKM value of WRKY70 was 281 at stage G, and was down-regulated to 130 at stage R and to 0.65 at stage C. Among the DE-lncRNAs that were co-expressed with WRKY70, only LNC_003467 were negatively correlated, LNC_001263 was highly expressed at stage G and was significantly down-regulated at stage C (Fig. 6). We speculated that WRKY70 might be a negative regulatory factor in the climacteric process, and that it might interact with LNC_001263. LNC_000987, LNC_001323, and LNC_003610 were up-regulated and highly expressed at stage C, and were co-expressed with TF3A, which participates in nucleic acid binding. As an important plant-specific TF, the AP2/ERF gene family plays crucial roles in plant growth, development, and stress responses. ERF genes mediate transcription of ethylene-regulated genes. The AP2/ERF family has been widely studied in many plant species, such as Arabidopsis thaliana  and Vitis vinifera . Among the DE-ERFs that we detected, ERF12, ERF76, and ERF80 were highly expressed at stages G and R, and were significantly down-regulated at stage C. Thus, ethylene-responsive genes may not be induced by these ERF genes and a portion of the ethylene synthesized may not be consumed in stage C. The concurrent significant up-regulation of CmACO1 and CmACO7 would notably promote ethylene biosynthesis. These findings suggested that the climacteric in C. melo might depend on the synthesis of a large amount of ethylene as well as reduction in consumption of ethylene.
The climacteric is a crucial physiological process during ripening of C. melo fruit. By performing high-throughput RNA sequencing and genome mapping, we constructed an informative dataset of lncRNAs and differentially expressed mRNAs associated with C. melo fruit ripening. The DE-lncRNAs and DE-mRNAs involved in plant hormone signal transduction, ethylene biosynthesis, sucrose biosynthesis and several transcription factors were investigated. Among these processes and metabolic pathways, the expression patterns of genes (CmACO, CmSAPK, CmETR), lncRNAs (LNC_000987, LNC_000693, LNC_001323, LNC_003610 and LNC_003380) and transcription factors (ERFs, NAC56, NAC73, and WRKY70) suggest that they may be closely associated with fruit ripening and the climacteric in C. melo. The biosynthesis and consumption of ethylene is likely to be a key element in the ripening process. The present results provide insight into the molecular mechanism of the respiratory climacteric, and lay a foundation for future studies of the ripening process in climacteric fruit of C. melo.
Plant material, RNA extraction and transcriptome sequencing
C. melo (Cultivar. Hetao) plants were grown in open experimental field (Bayan Nur, Inner-mongolia) from May to August in 2016 (6 °C–34 °C), temperature difference between day and night were more than 10 °C, self-pollination strictly. To identify expressed lncRNAs of Cucumis melo fruit, samples were collected at four developmental stages, namely growing (G), ripening (R), climacteric (C), and post-climacteric (P), with three biological replicates per stage. The G and R stage samples were harvested at 18 and 36 DAA, respectively. The collection of samples from harvested C stage fruit coincided with the peak respiration rate at 20 and 22 h post-harvest. The P stage samples were collected from post-climacteric harvested fruit at 48 h post-harvest. The fruit pressure was tested using fruit pressure tester (EFFEGIDI, FT 011, Italy). Four symmetry points were measured at the vertical section of the mesocarp. And soluble solids content was tested strictly following the instructions (Pocket refractometer: ATAGO, Cat. No. 3830, Japan). Respiration rate was also measured using a Fruit and Vegetable Respiratory Meter (3051H, TOP Instrument Co., Zhejiang, China). The time course of changes in the respiration rate of fruit harvested at 36 DAA.
After harvesting, the pericarp was immediately dissected, the flesh was frozen in liquid nitrogen, and stored at − 80 °C. Total RNA was isolated using the TRIzon Reagent (Kangwei Biotech, Beijing, China) in accordance with the manufacturer’s instructions. Genomic DNA was removed from extracted total RNA by treatment with DNase I (RNase-free; Takara Bio, Shiga, Japan). The RNA integrity, quality, and quantity were checked after 1% agarose gel electrophoresis with a NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, USA) and using the Nano 6000 Assay Kit with the Agilent 2100 Bioanalyzer™ system (Agilent Technologies, Santa Clara, CA, USA), respectively. Finally, 12 libraries were constructed and sequenced on an Illumina HiSeq 2500 platform.
Data filtering and genome mapping
Clean reads were obtained by removing reads containing adapter and poly-N sequences. Data processing of raw reads was quality checked trimmed for low-quality bases and adaptors by using Illumiprocessor Version 2.0 (https://illumiprocessor.readthedocs.io/en/latest/). The Q20 and Q30 percentages and GC content of the clean data were calculated. All downstream analyses were based on the clean data consisting of high-quality reads. Reference genome and gene model annotation files were downloaded from the Melonomics website (Genome/Melon_genome_v3.5.1, http://melonomics.cragenomica.es/). An index of the reference genome was built using Bowtie (Version 2.0.6) and paired-end clean reads were aligned to the reference genome using TopHat (Version 2.0.9, N = 2, library type = fr-first strand).
lncRNAs library construction and gene expression quantification
Transcripts shorter than 200 bp in length were filtered out, number of exons ≥2, expression levels FPKM≥0.5 were retained. Then, the transcripts were identified as lncRNAs or mRNAs by using the software programs CPC , and the Pfam database , the RNAs that lacked coding potential were considered to be candidates for lncRNAs and were used in subsequent analyses. The basic features of the obtained lncRNAs were characterized by comparison with mRNAs.
Cuffdiff (Version 2.1.1) was used to calculate fragments per kilo-base of FPKMs of both lncRNAs and coding genes. The FPKMs were computed by summing the FPKMs of transcripts in each gene group. Transcription with P-just < 0.05 were classified as differentially expressed.
Target gene prediction and enrichment analysis
Potential target genes of the lncRNAs were predicted according to their regulatory patterns, which were divided into cis- and trans-acting groups. For cis-regulator prediction, a lncRNA that was 100 kb up- or down-stream of the coding gene was defined as a co-located lncRNA. Prediction of trans-regulators was determined by the expression level; the correlation in expression between lncRNAs and coding genes was evaluated using the Pearson’s correlation coefficient (r > 0.95 or < − 0.95 and p < 0.05). The intersection between the two groups were extracted (Additional file 9). To gain insight into the functions of the lncRNAs and their corresponding target genes, GO terms and KEGG pathway enrichment analyses of the lncRNAs were performed. The GO ontologies were assigned using Blast2GO , and GO terms with a corrected P-value < 0.05 (Correction method, FDR) were considered significantly enriched. KEGG annotation was performed using KOBAS (Version 2.0, http://www.genome.jp/kegg/).
Systematic cluster analysis, protein–protein interaction and co-expression network construction
Systematic cluster analysis of DE-mRNAs, DE-TFs, and their targeted DE-lncRNAs was performed. For DE-mRNAs, those that were correlated with response to hormone stimulus, carotenoid or ethylene biosynthesis, sugar and main organic acid metabolism were selected. Transcription factors that showed fold change > 2 and FPKM > 30 were selected for the analysis. A heat map was generated using Cluster software (Version 3.0) and visualized using Java Treeview (Version 1.0.4). Protein–PPI analysis was based on the STRING database (http://string-db.org/), and a network was constructed by extracting the target gene list from the database. Target genes were aligned to the selected reference protein sequences using Blastx (Version 2.2.28), and the network was constructed in accordance with the known interaction of the reference species (Cucumis sativus). The gene IDs of the selected genes (including DE-TFs) were used to isolate the co-expression relations among target genes and lncRNAs, and the P-P interactions in PPI analysis results, respectively. Then, organized in a spreadsheet (Additional file 10) and constructed a co-expression network using Cytoscape (Version 3.7.0).
Quantitative real-time PCR of lncRNAs
Ten lncRNAs with only one transcript were selected and quantitative real-time PCR (qRT-PCR) was performed as validation. RNA isolation was performed as already described. Poly (A) Polymerase (Takara Bio, Shiga, Japan) was used to add poly A tails following the manufacturer’s instructions, then reverse transcription was performed using the PrimeScript™ RT reagent Kit with gDNA Eraser (Takara Bio, Shiga, Japan). qRT-PCR was performed using SYBR® Premix Ex Taq™ II (Takara Bio, Shiga, Japan) on a BioRad CFX96 instrument with default parameters. Relative gene expression levels were normalized to GAPDH and calculated using the 2−ΔΔCt method. Primers used for the validation experiments are shown in Table 2.
Availability of data and materials
The data has been submitting to the SRA (https://www.ncbi.nlm.nih.gov/sra), SRA accession: PRJNA543288.
Coding Potential Calculator
Days after anthesis
Differentially expressed lncRNAs
Auxin indole-3-acetic acid
Sucrose synthases Y
Wilhelm BT, Marguerat S, Watt S, Schubert F, Wood V, Goodhead I, Penkett CJ, Rogers J, Bähler J. Dynamic repertoire of a eukaryotic transcriptome surveyed at single-nucleotide resolution. Nature. 2008;453:1239.
The EPC, Birney E, Stamatoyannopoulos JA, Dutta A, Guigó R, Gingeras TR, Margulies EH, Weng Z, Snyder M, Dermitzakis ET, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature. 2007;447:799.
Igor U, Bartel DP. lincRNAs: genomics, evolution, and mechanisms. Cell. 2013;154(1):26–46.
Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP. Chromatin signature reveals over a thousand highly conserved large non-coding rnas in mammals. Nature. 2009;458(7235):223.
Jun L, Choonkyun J, Jun X, Huan W, Shulin D, Lucia B, Catalina A-H, Nam-Hai C. Genome-wide analysis uncovers regulation of long intergenic noncoding RNAs in Arabidopsis. Plant Cell. 2012;24(11):4333–45.
Wierzbicki AT, Ream TS, Haag JR, Pikaard CS. RNA polymerase V transcription guides ARGONAUTE4 to chromatin. Nat Genet. 2009;41(5):630–4.
Pruneski JA, Hainer SJ, Petrov KO, Martens JA. The Paf1 Complex Represses SER3 Transcription in Saccharomyces cerevisiae by facilitating intergenic transcription-dependent nucleosome occupancy of the SER3 promoter. Eukaryotic Cell. 2011;10(10):1283.
Wilusz JE. Long noncoding RNAs: re-writing dogmas of RNA processing and stability. Biochim Biophys Acta. 2016;1859(1):128–38.
Hadjiargyrou M, Delihas N. The intertwining of transposable elements and non-coding RNAs. Int J Mol Sci. 2013;14(7):13307–28.
Au PCK, Dennis ES, Wang M-B. Analysis of Argonaute 4-associated long non-coding RNA in Arabidopsis thaliana sheds novel insights into gene regulation through RNA-directed DNA methylation. Genes. 2017;8(8):198.
Ard R, Tong P, Allshire RC. Long non-coding RNA-mediated transcriptional interference of a permease gene confers drug tolerance in fission yeast. Nat Commun. 2014;5:5576.
Jain P, Sharma V, Dubey H, Singh PK, Kapoor R, Kumari M, Singh J, Pawar DV, Bisht D, Solanke AU, et al. Identification of long non-coding RNA in rice lines resistant to Rice blast pathogen Maganaporthe oryzae. Bioinformation. 2017;13(8):249–55.
Boerner S, McGinnis KM. Computational identification and functional predictions of long noncoding RNA in Zea mays. PLoS One. 2012;7(8):e43047.
Zhu B, Yang Y, Li R, Fu D, Wen L, Luo Y, Zhu H. RNA sequencing and functional analysis implicate the regulatory role of long non-coding RNAs in tomato fruit ripening. J Exp Bot. 2015;66(15):4483–95.
Deng F, Zhang X, Wang W, Yuan R, Shen F. Identification of Gossypium hirsutum long non-coding RNAs (lncRNAs) under salt stress. BMC Plant Biol. 2018;18(1):23.
Swiezewski S, Liu F, Magusin A, Dean C. Cold-induced silencing by long antisense transcripts of an Arabidopsis Polycomb target. Nature. 2009;462:799.
Cui J, Jiang N, Meng J, Yang G, Liu W, Zhou X, Ma N, Hou X, Luan Y. LncRNA33732-respiratory burst oxidase module associated with WRKY1 in tomato- Phytophthora infestans interactions. Plant J. 2018;97(5).
Saladié M, Cañizares J, Phillips MA, Rodriguez-Concepcion M, Larrigaudière C, Gibon Y, Stitt M, Lunn JE, Garcia-Mas J. Comparative transcriptional profiling analysis of developing melon (Cucumis melo L.) fruit from climacteric and non-climacteric varieties. BMC Genomics. 2015;16(1):440.
Boualem A, Fergany M, Fernandez R, Troadec C, Martin A, Morin H, Sari M-A, Collin F, Flowers JM, Pitrat M, et al. A conserved mutation in an ethylene biosynthesis enzyme leads to andromonoecy in melons. Science. 2008;321(5890):836–8.
Martin A, Troadec C, Boualem A, Rajab M, Fernandez R, Morin H, Pitrat M, Dogimont C, Bendahmane A. A transposon-induced epigenetic change leads to sex determination in melon. Nature. 2009;461:1135.
Giovannoni JJ. Fruit ripening mutants yield insights into ripening control. Curr Opin Plant Biol. 2007;10(3):283–9.
Pech JC, Bouzayen M, Latché A. Climacteric fruit ripening: Ethylene-dependent and independent regulation of ripening pathways in melon fruit. Plant Science. 2008;175(1):114–20.
Ayub R, Guis M, Amor MB, Gillot L, Roustan J-P, Latché A, Bouzayen M, Pech J-C. Expression of ACC oxidase antisense gene inhibits ripening of cantaloupe melon fruits. Nat Biotechnol. 1996;14(7):862–6.
Osorio S, Scossa F, Fernie AR. Molecular regulation of fruit ripening. Front Plant Sci. 2013;4:198.
Christophe P, Maricarmen GJ, Lynda H, Catherine D, Jean-Claude P, Alain L, Michel P, et al. Molecular and genetic characterization of a non-climacteric phenotype in melon reveals two loci conferring altered ethylene response in fruit. Plant Physiology. 2002;129(1):300–9.
Guo X, Xu J, Cui X, Chen H, Qi H. iTRAQ-based protein profiling and fruit quality changes at different development stages of oriental melon. BMC Plant Biol. 2017;17(1):28.
Zhang H, Wang H, Yi H, Zhai W, Wang G, Fu Q. Transcriptome profiling of Cucumis melo fruit development and ripening. Horticulture research. 2016;3:–16014.
Hao J, Niu Y, Yang B, Gao F, Zhang L, Wang J, Hasi A. Transformation of a marker-free and vector-free antisense ACC oxidase gene cassette into melon via the pollen-tube pathway. Biotechnol Lett. 2011;33(1):55–61.
Liu X, Hao L, Li D, Zhu L, Hu S. Long non-coding RNAs and their biological roles in plants. Genomics, proteomics & bioinformatics. 2015;13(3):137–47.
Winkel-Shirley B. It takes a garden. How work on diverse plant species has contributed to an understanding of flavonoid metabolism. Plant Physiol. 2001;127(4):1399–404.
Carrari F, Fernie AR. Metabolic regulation underlying tomato fruit development. J Exp Bot. 2006;57(9):1883–97.
Pattison RJ, Csukasi F, Catala C. Mechanisms regulating auxin action during fruit development. Physiol Plant. 2014;151(1):62–72.
Jia H, Jiu S, Zhang C, Wang C, Tariq P, Liu Z, Wang B, Cui L, Fang J. Abscisic acid and sucrose regulate tomato and strawberry fruit ripening through the abscisic acid-stress-ripening transcription factor. Plant Biotechnol J. 2016;14(10):2045–65.
Chang C, Kwok SF, Bleecker AB, Meyerowitz EM. Arabidopsis ethylene-response gene ETR1: similarity of product to two-component regulators. Science. 1993;262(5133):539.
Abdullah M, Cao Y, Cheng X, Meng D, Chen Y, Shakoor A, Gao J, Cai Y. The sucrose synthase gene family in Chinese pear (Pyrus bretschneideri Rehd.): structure, expression, and evolution. Molecules (Basel, Switzerland). 2018;23(5):1144.
Liu M, Pirrello J, Chervin C, Roustan J-P, Bouzayen M. Ethylene control of fruit ripening: revisiting the complex network of transcriptional regulation. Plant Physiol. 2015;169(4):2380–90.
Costa F, Alba R, Schouten H, Soglio V, Gianfranceschi L, Serra S, Musacchi S, Sansavini S, Costa G, Fei Z, et al. Use of homologous and heterologous gene expression profiling tools to characterize transcription dynamics during apple fruit maturation and ripening. BMC Plant Biol. 2010;10:229.
Lü P, Yu S, Zhu N, Chen Y-R, Zhou B, Pan Y, Tzeng D, Fabi JP, Argyris J, Garcia-Mas J, et al. Genome encode analyses reveal the basis of convergent evolution of fleshy fruit ripening. Nature Plants. 2018;4(10):784–91.
Zhu M, Chen G, Zhou S, Tu Y, Wang Y, Dong T, Hu Z. A new tomato NAC (NAM/ATAF1/2/CUC2) transcription factor, SlNAC4, functions as a positive regulator of fruit ripening and carotenoid accumulation. Plant Cell Physiol. 2014;55(1):119–35.
Rios P, Argyris J, Vegas J, Leida C, Kenigswald M, Tzuri G, Troadec C, Bendahmane A, Katzir N, Pico B, et al. ETHQV6.3 is involved in melon climacteric fruit ripening and is encoded by a NAC domain transcription factor. Plant J. 2017;91(4):671–83.
An JP, Zhang XW, Bi SQ, You CX, Wang XF, Hao YJ. MdbHLH93, an apple activator regulating leaf senescence, is regulated by ABA and MdBT2 in antagonistic ways. New Phytol. 2019;222(2).
Sun C, Palmqvist S, Olsson H, Borén M, Ahlandsberg S, Jansson C. A novel WRKY transcription factor, SUSIBA2, participates in sugar signaling in barley by binding to the sugar-responsive elements of the iso1 promoter. Plant Cell. 2003;15(9):2076–92.
Zhou H, Li Y, Zhang Q, Ren S, Shen Y, Qin L, Xing Y. Genome-wide analysis of the expression of WRKY family genes in different developmental stages of wild strawberry (Fragaria vesca) fruit. PLoS One. 2016;11(5):e0154312.
Nakano T, Suzuki K, Fujimura T, Shinshi H. Genome-wide analysis of the ERF gene family in Arabidopsis and rice. Plant Physiol. 2006;140(2):411–32.
Licausi F, Giorgi FM, Zenoni S, Osti F, Pezzotti M, Perata P. Genomic and transcriptomic analysis of the AP2/ERF superfamily in Vitis vinifera. BMC Genomics. 2010;11:719.
Kong L, Zhang Y, Ye Z-Q, Liu X-Q, Zhao S-Q, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research. 2007;35(suppl_2):W345–9.
Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, Potter SC, Punta M, Qureshi M, Sangrador-Vegas A, et al. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44(D1):D279–85.
Conesa A, Terol J, García-Gómez JM, Talón M, Robles M, Götz S. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6 %@ 1367–4803.
We thank Robert McKenzie, PhD, from Liwen Bianji, Edanz Group China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.
This work was supported by the National Natural Science Foundation of China (No. 31560561 and 31360486) and the Inner Mongolia Natural Science Foundation of China (No. 2017ZD05). The funding bodies were not involved in the design of the study, data collection, interpretation of data, or in writing the manuscript.
Ethics approval and consent to participate
These plant materials are widely used all over the world and no permits are required for the collection of plant samples. The seeds are preserved by the supervisor of our lab for decades. The plant materials are maintained in accordance with the institutional guidelines of the College of Life Sciences, Inner Mongolia University, China. This article did not contain any studies with human participants or animals and did not involve any endangered or protected species.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Details of lncRNAs. (XLSX 275 kb)
Cis- and Trans-interaction between lncRNAs and mRNAs. (XLSX 20120 kb)
The GO terms of co-located target genes. (XLS 34 kb)
Top 20 enriched KEGG pathways of co-located target genes. (XLS 33 kb)
Enriched cis-target genes in plant hormone signal transduction and Carotenoid biosynthesis. (XLS 191 kb)
GO terms of co-expressed genes. (XLS 2482 kb)
Top 20 enriched KEGG pathways of co-expressed target genes. (XLS 63 kb)
Gene ID of the genes were discussed. (XLSX 31 kb)
The intersection between co-location and co-expression groups. (XLSX 60 kb)
Details of the co-expression networks. (XLSX 551 kb)
Figure S1. Clustering analysis of all the differentially expressed lncRNAs. (PDF 149 kb)
Figure S2. The results of the RT-qPCR confirmed the expression patterns of the selected lncRNAs were consistent with the expression levels calculated from the RNA-seq data. (JPG 200 kb)
About this article
Cite this article
Tian, Y., Bai, S., Dang, Z. et al. Genome-wide identification and characterization of long non-coding RNAs involved in fruit ripening and the climacteric in Cucumis melo. BMC Plant Biol 19, 369 (2019). https://doi.org/10.1186/s12870-019-1942-4