Melonet-DB, a Grand RNA-Seq Gene Expression Atlas in Melon (Cucumis melo L.)

Abstract

Melon (Cucumis melo L.) is an important Cucurbitaceae crop produced worldwide, exhibiting wide genetic variations and comprising both climacteric and non-climacteric fruit types. The muskmelon cultivar “‘Earl’s favorite Harukei-3 (Harukei-3)”’ known for its sweetness and rich aroma is used for breeding of high-grade muskmelon in Japan. We conducted RNA sequencing (RNA-seq) transcriptome studies in 30 different tissues of the ‘Harukei-3’ melon. These included root, stems, leaves, flowers, regenerating callus and ovaries, in addition to the flesh and peel sampled at seven stages of fruit development. The expression patterns of 20,752 genes were determined with fragments per kilobase of transcript per million fragments sequenced (FPKM) >1 in at least one tissue. Principal component analysis distinguished 30 melon tissues based on the global gene expression profile and, further, the weighted gene correlation network analysis classified melon genes into 45 distinct coexpression groups. Some coexpression groups exhibited tissue-specific gene expression. Furthermore, we developed and published web application tools designated “‘Gene expression map viewer”’ and “‘Coexpression viewer”’ on our website Melonet-DB (http://melonet-db.agbi.tsukuba.ac.jp/) to promote functional genomics research in melon. By using both tools, we analyzed melon homologs of tomato fruit ripening regulators such as E8, RIPENING-INHIBITOR (RIN) and NON-RIPENING (NOR). The “‘Coexpression viewer”’ clearly distinguished fruit ripening-associated melon RIN/NOR/CNR homologs from those expressed in other tissues. In addition, several other MADS-box, NAM/ATAF/CUC (NAC) and homeobox transcription factor genes were identified as fruit ripening-associated genes. Our tools provide useful information for research not only on melon but also on other fleshy fruit plants.

Introduction

With rapidly expanding knowledge of gene function such as the mutant phenotypes, gene expression profiles and protein structures, bioinformatics databases have become integral to routine research in genomics, genetics and molecular biology in both model and non-model biological entities. Due to the cost reductions achieved in microarray and next-generation sequencing (NGS) technologies, genome-wide omics datasets are being generated at higher frequencies than before. These emerging changes pose a huge challenge in developing comprehensive information databases and implementing user-friendly interactive web tools to enhance our understanding of the biological systems.

During the 1990s, with the introduction of genome science, whole-genome information databases became available, developed mainly from the model biological organisms such as Arabidopsis. Since the Arabidopsis genome was the first published plant genome, its gene database has evolved to be highly comprehensive. Although the Arabidopsis Information Resource (TAIR; https://www.arabidopsis.org/) plays a central role in providing information related to DNA/protein sequence, mutant phenotype and relevant literature, some databases have been developed and specialized for the analysis of transcriptome information. For instance, the Arabidopsis eFP browser was the first web tool that enables in silico gene expression analysis (Winter et al. 2007). A number of transcriptome datasets have been integrated in the eFP browser, including those of plant tissues and organs, seed development and germination, abiotic/biotic stress conditions, phytohormone responses and natural variations. The expansion of such transcriptome datasets has also enabled genome-wide gene expression correlation network analysis (Obayashi et al. 2007, Obayashi et al. 2009, Usadel et al. 2009). Investigation of gene to gene correlation at the whole-genome level facilitates determination of coexpression groups (CGs) or genes that are involved in a specific biological process (Fukushima et al. 2008, Bassel et al. 2011, Lin et al. 2011, Mantegazza et al. 2014, Nagel et al. 2015, Silva et al. 2016). In particular, coexpression studies are powerful when combined with metabolite datasets; both transcription factor and enzyme genes have been identified through this strategy (Hirai et al. 2007, Yonekura-Sakakibara et al. 2007, Yonekura-Sakakibara et al. 2008, Sakurai et al. 2011). Transcriptome studies have also been widely conducted in non-model plants, including economically important crops such as rice, soybean and tomato (He et al. 2009, Guttikonda et al. 2010, Ozaki et al. 2010, Cao et al. 2012, Fukushima et al. 2012, Joshi et al. 2012, Sato et al. 2013, Yu et al. 2014). In tomato, the model of fleshy plants, “‘TOMATOMICS”’ and “‘Sol Genomics Network”’, have been developed as comprehensive tomato genome databases (Fernandez-Pozo et al. 2015, Kudo et al. 2017). These enable digital gene expression analysis as well as correlation network analysis. The ATTED II database contains coexpression datasets of nine plant species, including Arabidopsis and tomato (Aoki et al. 2016). Similar to Arabidopsis, coexpression studies appear to be effective in identifying genes involved in secondary metabolism; in Medicago truncatula, a cytochrome P450 gene that is involved in triterpenoid biosynthesis has been identified through a survey of the coexpression database (Fukushima et al. 2011). In Solanaceae crops, coexpression studies combining transcriptome and metabolite datasets identified key genes involved in toxic steroidal glycoalkaloids such as α-solanine (Itkin et al. 2013).

Melon (Cucumis melo L.) is an economically important fruit crop of the Cucurbitaceae family that is produced worldwide. Although the melon fruit is an important source of vitamins and minerals for humans, it exhibits a wide range of natural variation, including climacteric (ethylene-producing; var. cantalupensis and reticulatus) and non-climacteric (non-ethylene-producing; var. inodorus) fruit ripening types (Kendall and Ng 1988, Stepansky et al. 1999, Moreno et al. 2008, Diaz et al. 2011, Vegas et al. 2013, Monforte et al. 2014). This fruit also exhibits great variation in fruit shape, sugar content and disease resistance ability (Monforte et al. 2004, Fukino et al. 2008, Obando-Ulloa et al. 2009, Harel-Beja et al. 2010, Díaz et al. 2014). The first whole-genome information on melon was published based on the homozygous double-haploid line ‘DHL92’, which was derived from a cross between the ‘Piel de Sapo’ (var. inodorus) and the Korean melon ‘Songwhan Charmi’ (Garcia-Mas et al. 2012). A number of transcriptome studies have been conducted on several distinct melon cultivars from different geographic origins (Saladie et al. 2015, Li et al. 2016, Zhang et al. 2016, Shin et al. 2017). A microarray study comparing global gene expression profiles of climacteric and non-climacteric melon cultivars has identified potential candidate genes that may account for differences observed in the fruit ripening types (Saladie et al. 2015). However, as most transcriptome studies on melon have focused on fruit ripening, and did not include non-fruit organs or tissues, the nature of fruit ripening-associated genes is not fully understood.

In this study, we present a large gene expression atlas of melon that was constructed from RNA sequencing (RNA-seq) data obtained from 30 different tissues, including fruit developmental stages from days after flowering (DAF) 0–50. For this purpose, we used Japanese muskmelon ‘Earls favorite Harukei-3’ (‘Harukei-3’) because this melon is important for breeding of expensive muskmelon in Japan. Based on the global gene expression profiles of 30 tissues, we could successfully classify melon genes into different CGs. This attempt identified a number of fruit ripening-associated genes. To enable in silico gene expression analysis in melon, we further constructed web application tools designated “‘Gene expression map viewer”’ and “‘Coexpression viewer”’ and posted them on our website Melonet-DB (http://melonet-db.agbi.tsukuba.ac.jp/). The former tool enables digital gene expression comparison from multiple queries, while the latter tool performs coexpression network analysis. Here, we first report the result of the RNA-seq study [i.e. weighted gene correlation network analysis (WGCNA) clustering], and then demonstrate the usability of our tools through a case study of fruit ripening-associated genes in melon. Although the fruit ripening melon quantitative trait locus (QTL) ETHQV6.3 (MELO3C016540) was recently shown to encode a NAC domain transcription factor that is similar to that ecoded by the tomato NON-RIPENING (NOR) gene (Rios et al. 2017), the entire gene regulatory network associated with melon fruit ripening remains unclear. We show that our database tools identified several candidate melon homologs for tomato fruit ripening regulators, including RIPENING-INHIBITOR (RIN), NOR and COLORLESS NON-RIPENING (CNR).

Results and Discussion

Global gene expression profiles of 30 different tissues of the ‘Harukei-3’ melon

In order to obtain global gene expression profile information, the ‘Harukei-3’ melon plants were grown under greenhouse conditions from early spring to summer and then RNA-seq analysis was conducted in 30 different tissues. These tissues include fully expanded leaves (sixth, ninth and 12th true leaves), young leaves (<5 cm length), stems (middle and upper part), root, tendrils, petals, anthers, female flowers (stigma, ovary at DAF 0, 2 and 4), fruits (flesh and epicarp, DAF 8, 15, 22, 29, 36, 43 and 50), dry seeds and calli (Fig. 1). Except for fruit tissues, all tissues were obtained just after male flower formation. In general, the timing of fruit harvest varies between melon cultivars dependent on fruit ripening behavior. While climacteric melon fruits are harvested around 35 DAF (i.e. ‘Vedrantais’), the harvest time of non-climacteric melon is around 70 DAF (i.e. ‘Honeydew’). Fruits of the ‘Harukei-3’ melon are generally harvested around DAF 50. Under our experimental conditions, ethylene emission was not observed at DAF 50. However, >60% of fruits exhibited ethylene emission within 3 weeks of subsequent post-harvest storage (Supplementary Fig. S1), indicating that ethylene emission occurs at a much later time in ‘Harukei-3’ when compared with typical climacteric melon such as ‘Vedrantais’.

Fig. 1

Melon plant tissue samples used for RNA-seq study. (A) A cartoon illustrating 30 ‘Harukei-3’ melon plant tissue samples used for RNA-seq analysis. These include true leaves (sixth, ninth and 12th), young leaves (<5 cm length), stems (middle and upper part), root, tendril, petal, anther, female flower (stigma and ovary), fruit (flesh and epicarp, DAF 2, 4, 8, 15, 22, 29, 36, 43 and 50), dry seed and callus. All tissues were obtained at the timing of male flower formation except fruit tissues. (B) Melon plants and fruits grown under greenhouse conditions. Photographs of melon fruits were taken at DAF 2, 4, 8, 15, 22, 29, 36, 43 and 50 (right). Scale bar = 1 cm.

Using the Illumina Hiseq NGS platform, we obtained a total of 34 Gb of paired-end short read data (after quality control and removal of microbial contamination reads). The CM3.5.1 melon genome information (Garcia-Mas et al. 2012) was used as a reference, and NGS reads were aligned to the reference at an average rate of 91.6% (Fig. 2A). Of the 27,427 annotated genes, 20,752 genes exhibited fragments per kilobase of transcript per million fragments sequenced (FPKM) > 1 for at least one tissue. When the expression dataset of 20,752 genes was subjected to principal component analysis (PCA), 30 melon tissues were distinguished on the PCA plot (Fig. 2B). In particular, gene expression profiles were clearly differentiated between fruit tissues, leaves and other sink tissues such as dry seed and anther. Early fruit developmental stages (DAF 15–29) and late fruit developmental stages (DAF 43 and 50) were also distinguished, suggesting differences in gene expression profiles between these stages.

Fig. 2

Principal component analysis (PCA) of the global gene expression profile in 30 ‘Harukei-3’ melon tissues. (A) RNA-seq read alignment against the CM3.5.1 melon genome reference. Frequency of aligned reads (top) or the number of genes with FPKM >1 (bottom) are shown. (B) PCA based on the global gene expression profiles of 30 melon tissues. Log₂-transformed values were subjected to PCA implemented in R software. Comp 1 and 2 explain 70.1% of total variance. Ovaries (DAF 0, 2 and 4), expanded leaves (sixth, ninth and 12th true leaves), and fruit tissues at early development stages (DAF 15, 22 and 29) or late development stages (DAF 43 and 50) are indicated by blue, green, orange or red circles, respectively. Ant, anther; L6, sixth true leaves; L9, ninth true leaves; L12, 12th true leaves; CL, regenerating callus; YL, young leaves; US, upper side stem (elongating zone); MS, middle zone stem; Tn, tendril; Pet, flower petal; Ov0, female flower ovary without pollination; Ov2, female flower ovary at DAF 2; Ov4, female flower ovary at DAF 4; DS, dry seed; Rt, root; FL, fruit flesh at certain DAFs; EP, fruit epicarp at certain DAFs.

Identification of coexpression groups by genome-wide correlation network analysis

Next, we performed WGCNA to identify CGs in the global gene expression dataset (Langfelder and Horvath 2008). According to the results, 17,597 genes were classified into 45 CGs when the P-value threshold was 1^e–3 (Fig. 3A;Table 1). Some CGs exhibited clear tissue specificity. For instance, gene expression of CG15, CG16, CG17, CG19, CG27, CG29, CG30, CG23, CG24, CG32 and CG33 was specific to, or was the strongest in dry seed, stigma, petal, anther, root, callus, tendril, young leaves, expanded leaves, the middle part of the stem and the upper part of the stem, respectively (Fig. 3B). Strong expression of CG36, CG38 and CG37 was observed in the ovary of female flowers at DAF 0, 2 and 4, respectively, suggesting that genes classified into these CGs might have a functional role during the initial stages of fruit development. In addition, CG10, CG11, CG12, CG13 and CG14 exhibited fruit-specific gene expression, as their levels were found to be high from DAF 22 to DAF 50 in fruit flesh tissues, implying their role in fruit development and ripening. Phylogenetic tree analysis using eigengene information (artificial numeric vector that represents the expression profile of each CG) showed that CGs with similar expression patterns were clustered in the neighboring nodes (Fig. 3C). Some CGs showed strong gene expression not only in specific tissues but also in several distinct tissues. For example, CG35 and CG36 exhibited preferential gene expression in both stem and root. In addition, the expression of CG18 was strong in both anther and petal tissue samples, while that of CG29 was stronger in both the root and callus. It is likely that genes classified into such ‘tissue-overlapping’ CGs are involved in common molecular process(es) in different tissues. In addition to the set of 30 tissues, we also performed WGCNA in 29 tissues without regenerating callus sample because callus is not a naturally occurring plant tissue (Supplementary Fig. S2; Supplementary Table S1). In total, 20,650 genes had FPKM >1, and 40 CGs were identified by WGCNA. Similar to WGCNA in 30 tissues, several CGs exhibited tissue-specific gene expression patterns (Supplementary Fig. S2B, C).

Table 1

Co-expression group identification by WGCNA in 30 ‘Harukei-3’ melon tissues

Group	P < 1e-3	P < 1e-5	P < 1e-10	Tissue specificity
CG1	101	60	13	–
CG2	51	48	33	–
CG3	45	29	12	–
CG4	75	53	18	–
CG5	42	33	19	–
CG6	98	66	13	–
CG7	52	35	12	–
CG8	26	21	11	–
CG9	21	19	10	–
CG10	180	120	34	FL43, FL50
CG11	55	43	26	FL50
CG12	336	175	12	FL22
CG13	158	118	41	FL43
CG14	199	148	39	FL29, FL36, FL43
CG15	755	578	372	DS
CG16	222	164	76	Sti
CG17	578	455	203	Pet
CG18	284	232	116	Ant, Pet
CG19	1759	1421	855	Ant
CG20	26	21	12	–
CG21	117	91	47	–
CG22	44	42	25	–
CG23	295	271	174	YL
CG24	1933	1497	753	L6, L9, L12
CG25	605	458	229	L6, L9, L12, YL
CG26	363	327	173	YL, L6, L9, L12
CG27	764	596	362	Rt
CG28	782	615	363	CL
CG29	40	35	25	CL, Rt
CG30	165	142	83	Tn
CG31	66	63	41	Tn, Rt
CG32	106	86	40	MS
CG33	194	174	102	US
CG34	194	161	74	US, MS, Rt
CG35	87	78	46	Rt, US, MS
CG36	37	30	22	Ov0
CG37	49	41	22	Ov4
CG38	167	135	56	Ov2
CG39	2254	1475	304	–
CG40	2034	1446	297	–
CG41	1690	1293	436	Ov0-4, CL, YL, US
CG42	67	47	18	–
CG43	80	55	18	–
CG44	146	119	53	–
CG45	255	199	56	–

The probability of group membership for each gene is indicated by P-values that were calculated by WGCNA. In total, 17,597 genes were classified into 45 co-expression modules (P-value threshold = 1^e–3). Some co-expression groups exhibited tissue specificity as shown in Fig. 3B.

Ant, anther; L6, sixth true leaves; L9, ninth true leaves; L12, 12th true leaves; CL, regenerating callus; YL, young leaves; US, upper side stem (elongating zone); MS, middle zone stem; Tn, tendril; Pet, flower petal; Ov0, female flower ovary without pollination; Ov2, female flower ovary at DAF 2; Ov4, female flower ovary at DAF 4; DS, dry seed; Rt, root; FL, fruit flesh at certain DAFs; EP, fruit epicarp at certain DAFs.

Fig. 3

Coexpression group identification by WGCNA in 30 ‘Harukei-3’ melon tissues. (A) A heat map showing WGCNA clustering of melon genes. (B) Gene expression patterns of tissue-specific coexpression modules. P-values indicate the probability of group membership for each gene, and the cut-off threshold was set to 1^e–3. Ant, anther; L6, sixth true leaves; L9, ninth true leaves; L12, 12th true leaves; CL, regenerating callus; YL, young leaves; US, upper part of stem (elongating zone); MS, middle part of stem; Tn, tendril; Pet, flower petal; Ov0, female flower ovary without pollination; Ov2, female flower ovary at DAF 2; Ov4, female flower ovary at DAF 4; DS, dry seed; Rt, root; FL, fruit flesh at certain DAFs; EP, fruit epicarp at certain DAFs. (C) A phylogenetic tree showing the eigengene relationship. Eigengenes are artificial expression vectors that indicate representative expression profiles of each coexpression group. Coexpression modules CG10–CG14, CG15, CG16, CG17, CG18, CG19, CG23–CG26, CG27, CG28, CG29, CG30, CG31, CG32, CG33, CG34–G35 and CG36–CG38 exhibited the highest gene expression levels in fruit flesh, dry seed, stigma, petal, petal/anther, anther, expanded/young leaves, root, callus, callus/root, tendril, root/tendril, middle part of stem, upper part of stem, stems/root and ovary, respectively.

In order to evaluate the consistency of WGCNA clustering, we investigated keyword enrichment analysis using ID information from InterProScan (Jones et al. 2014). Of the 34,848 transcripts derived from 27,427 genes in the CM3.5.1 melon reference, 23,340 transcripts were found to encode proteins with an InterProScan ID (IPR ID) (Supplementary Data set S1). We focused on 11 CG categories that exhibited clear tissue-specific gene expression (i.e. anther, callus, dry seed, petal, root, stigma, tendril, fruit flesh, leaves, ovary and stem/root) (Supplementary Dataset S2), and found that some IPR IDs were over-represented in tissue-specific CGs (Fig. 4). For instance, pectin metabolism-related genes such as pectin lyase (IPR012334), pectin esterase inhibitor (IPR006501) and pectin esterase (IPR000070) were over-represented in anther-specific CGs (Fig. 4A; CG19 in Fig. 3B), suggesting that pectin metabolism-related genes are important for anther function (i.e. anther dehiscence and/or pollen dispersal) (Francis et al. 2006, Tian et al. 2006). Heat shock chaperon genes (IPR008978, IPR002068 and IPR013126) and seed storage protein genes (IPR000136, IPR007011, IPR005513, IPR018930 and IPR004238) were predominantly expressed in dry seed (Fig. 4B), consistent with their roles in protein stabilization, dehydration tolerance and/or nutrient storage in desiccated seed (Nakabayashi et al. 2005, Bassel et al. 2011). In root-specific CGs, lipoxygenase (IPR013819, IPR000907 and IPR027433) and disease resistance genes (IPR002182 and IPR004265) were over-represented compared with other tissues (Fig. 4C), suggesting that genes involved in microbial interaction and/or defense response are necessary in this tissue. Photosystem- (IPR023329, IPR022796 and IPR001344) and redox-related genes (IPR012336, IPR013766 and IPR005746) were also over-represented in leaf tissue samples. BURP protein domain-containing polygalacturonases were over-represented in fruit flesh tissues (Fig. 4E) (Park et al. 2015). Sieve element occlusion-related genes (IPR027944 and IPR027942) were found to be enriched in the CGs specific to stem and root (Fig. 4F), suggesting that they play an important role in translocation. Although the ‘Harukei-3’ plants used in this study were grown under fluctuating greenhouse conditions, these results suggest that melon genes have been successfully classified into CGs consistent with their gene functions.

Fig. 4

InterProScan ID enrichment analysis in tissue-specific coexpression groups. InterProScan ID was assigned to melon genes and ID enrichment was analyzed in tissue-specific coexpression modules. Enrichment of specific functional gene groups (i.e. pectin metabolizing genes, storage protein genes and photosynthetic genes) was observed for coexpression modules specific to anther (A), dry seed (B), root (C), leaves (D), fruit flesh (E) and stem/root (F).

Digital gene expression analysis by Melonet-DB ‘Gene expression map viewer’ and ‘Coexpression viewer’: a case study of fruit ripening-related genes

In order to simplify exploration of global gene expression information obtained in this study, we developed two web application tools designated “‘Gene expression map viewer”’ (http://melonet-db.agbi.tsukuba.ac.jp/cgi-bin/melview.cgi) and “‘Coexpression viewer”’ (http://melonet-db.agbi.tsukuba.ac.jp/cgi-bin/melcor.cgi) and published on our website “‘Melonet-DB”’. “‘Gene expression map viewer”’ enables digital gene expression comparison over 30 melon tissues through multiple queries (see instructions on our website). “‘Coexpression viewer”’ enables coexpression network analysis based on the WGCNA clustering results, and as described above.

Here, we present a case study of these web application tools. First, we analyzed the expression of volatile-related genes using the “‘Gene expression map viewer”’. Melon fruit is known for its excellent and rich aroma, and, particularly in muskmelon, aromatic volatile production is an important consumer preference, and is useful for marketing. In the melon genome, there are at least 12 alcohol dehydrogenase genes (CmADH1–CmADH12) and one ADH-like gene (CmFDH1) (Garcia-Mas et al. 2012, Chen et al. 2016, Jin et al. 2016). According to the ‘Harukei-3’ gene expression database, fruit ripening-specific gene expression was observed in CmADH1, CmADH2 and CmADH10 (Fig. 5). In contrast, the expression of other CmADH genes such as CmADH8 and CmADH12 was observed in all organs (Supplementary Fig. S3). In addition to ADH genes, alcohol acyl-transferase (AAT) genes have been also reported to be involved in aromatic volatile production in melon fruits (Shalit et al. 2001, Lucchetta et al. 2007). In the melon genome, there are at least four AAT genes (CmAAT1– CmAAT4). CmAAT1, which is involved in the production of ester compound such as cinnamyl acetate, has been reported to be associated with quantitative variation in aromatic volatiles in a recombinant inbred line population (El-Sharkawy et al. 2005, Freilich et al. 2015). Likewise, CmAAT1 exhibited a highly fruit ripening-specific gene expression pattern in the “‘Gene expression map viewer”’ (Fig. 5). Although the exact role of CmAAT2 is not fully understood, its gene expression was also found to be specific to fruit ripening. Comparison of digital gene expression images showed that the expression of CmADH1, CmADH2 and CmADH10 precedes that of CmAAT1 and CmADH2 (Fig. 5). As aroma production becomes evident at the later stages of fruit development, say from DAF 43 in ‘Harukei-3’, CmAAT1 and CmAAT2 are more likely to be linked to aroma production compared with CmADH genes. Unlike CmAAT1 and CmAAT2, the expression of CmAAT3 was not observed in fruit tissues, but was instead observed in other vegetative tissues, while CmAAT4 expression levels were below the detection limit in our RNA-seq study (Supplementary Fig. S3).

Fig. 5

Digital expression analysis of melon volatile-related genes by Melonet-DB’s ‘Gene expression map viewer’. Digital expression analysis of melon alcohol dehydrogenase (ADH) genes and alcohol acyl-transferase (AAT) genes was analyzed by Melonet-DB ‘Gene expression map viewer’ (http://melonet-db.agbi.tsukuba.ac.jp/cgi-bin/melview.cgi). Note the expression of MELO3C023685 (CmADH1), MELO3C014897 (CmADH2), MELO3C026554 (CmADH10), MELO3C024771 (CmAAT1) and MELO3C024766 (CmAAT2) is highly associated with fruit ripening in the ‘Harukei-3’ melon.

In order to identify candidate melon genes for fruit ripening regulators, we queried for melon homologs for tomato fruit ripening-related genes using the BLASTp tool. In tomato, fruit ripening is controlled by the concerted action of RIN, NOR, CNR and other transcription factor genes (Giovannoni 2004, Giovannoni 2007, Seymour et al. 2013, Karlova et al. 2014). RIN encodes a SEPALLATA-type MADS-box transcription factor, and several other MADS-box genes such as MADS1, TOMATO AGAMOUS-LIKE 1 (TAGL1) and FRUITFULL1/2 (FUL1/2) have been shown to be involved in tomato fruit ripening (Itkin et al. 2009, Vrebalov et al. 2009, Bemer et al. 2012, Dong et al. 2013, Shima et al. 2014). NOR encodes a NAC domain transcription factor (Tigchelaar et al. 1973, Giovannoni et al. 2004, Hurst and Loeffler 2014, Yuan et al. 2016), and another NAC transcription factor, NAC4, has been reported to play a role in fruit ripening (Zhu et al. 2014). CNR encodes a squamosa-promoter binding protein (SBP) domain transcription factor that was identified through epigenetic mutation (Eriksson et al. 2004). In addition to these, APETALA2a (AP2a) and HB-1 transcription factors have been identified and shown to be involved in the control of fruit ripening (Lin et al. 2008, Chung et al. 2010). Additionally, the tomato E8 gene is often used as a marker of fruit ripening. The E8 gene affects fruit ripening because co-suppression of the E8 gene results in increased ethylene production in the fruit (Kneissl and Deikman 1996). BLASTp analysis indicated that there are at least 38 non-redundant melon homologs for these tomato fruit ripening-related genes in the CM3.5.1 melon reference genome (Table 2). Among them, Rios et al. (2017) recently identified MELO3C016540 (CmNAC-NOR) as a fruit ripening QTL regulating climacteric ripening in melon. This gene encodes a NAC transcription factor similar to the tomato NOR. On the other hand, the role of the remaining melon homologs in fruit ripening is largely unknown (Table 2). All of these melon homologs had BLAST E-values >1^e–100 and protein sequence similarity <70%; therefore, it is difficult to determine the functional relationship between tomato and melon genes based merely on sequence similarity data. Thus, in order to limit fruit ripening-associated melon genes further, all of these melon homologs for tomato fruit ripening genes (E8, RIN, MADS1, TAGL1, FUL1/2, NOR, NAC4 and CNR) were subjected to coexpression network analysis using Melonet-DB “‘Coexpression viewer”’. Aromatic volatile-related genes (Fig. 5) were also combined in the query list as positive controls because it is difficult to find specific coexpression clusters without controls. As a result, 38 melon genes were classified into five clusters according to the coexpression patterns (Fig. 6A). Among these, of particular interest is the fruit-ripening associated coexpression cluster (Fig. 6A). This cluster contains not only CmADH1, CmADH2, CmADH10, CmADH1 and CmADH2 genes, but also two RIN/MADS1/TAGL1/FUL1/2 homologs (MELO3C026300 and MELO3C002691), two NOR homologs [MELO3C016540 (CmNAC-NOR) and MELO3C018237], one HB-1 like homolog and one AP2a-like homolog. Examination of the expression patterns of these transcription factor genes using “‘Gene expression map viewer”’ revealed their high expression levels in fruit flesh tissues, especially around DAF 29 (Fig. 7A). Interestingly, MELO3C026300, MELO3C002691, MELO3C016540 (CmNAC-NOR) and MELO3C007572 were the best hit homologs of tomato RIN/MADS1, TAGL1, NOR and AP2a, respectively (Table 2), suggesting key roles in melon fruit ripening. In addition to MELO3C016540 (CmNAC-NOR), MELO3C018237 also exhibited significant homology to tomato NOR (E-value = 6.41E-100), implying the existence of additional NOR homologs in the melon genome. Although the best hit homologs for tomato E8 and HB-1 were not classified in the fruit ripening-associated cluster (Fig. 6A), other homologs, MELO3C022978 (E8 homolog) and MELO3C021978 (HB-1 homolog), were identified as fruit ripening-associated genes. In addition to this cluster, another cluster appeared to have fruit development-associated genes (Fig. 6A). According to the “‘Gene expression map viewer”’, the expression of CNR homologs (MELO3C025597, MELO3C002618 and MELO3C009639) and TAGL1/FUL1/2 homologs (MELO3C022209 and MELO3C011409) was found to be strong during the early stages of fruit development from DAF 0 to 4 (Fig. 7B). Although there were no apparent melon homologs for tomato CNR (E-value >1^e–40; ∼60% protein sequence similarity), it is possible that either MELO3C025597, MELO3C002618 or MELO3C009639 may have functions similar to tomato CNR.

Table 2

Melon homologs for tomato fruit ripening-related genes

Tomato gene (query)	Name	Type	Hit melon gene	Name	E-value	Score	%identity
Solyc09g089580.2.1	E8	Fe(II)/2OG-dependent oxygenase	MELO3C006439	E8-L1	3.67E-131	379	52.1
			MELO3C006438	E8-L2	9.51E-129	374	50.3
			MELO3C022978	E8-L3	1.26E-112	332	43.1
			MELO3C026436	E8-L4	4.78E-89	273	42.4
			MELO3C016225	E8-L5	4.11E-58	192	35.3
Solyc03g044300.2.1	AP2a	AP2	MELO3C007572	AP2a-L1	2.78E-121	355	62.1
			MELO3C020848	AP2a-L2	7.99E-120	356	63.5
			MELO3C014261	AP2a-L3	1.52E-107	318	61.5
			MELO3C025726	AP2a-L4	5.32E-106	321	57.6
			MELO3C009307	AP2a-L5	4.99E-84	260	76.8
Solyc02g086930.2.1	HB1	Homeobox	MELO3C021534	HB-L1	1.91E-95	282	58.0
			MELO3C009948	HB-L2	1.32E-91	272	56.3
			MELO3C002209	HB-L3	1.28E-43	151	60.4
			MELO3C021978	HB-L4	1.79E-40	142	62.4
			MELO3C007736	HB-L5	8.06E-40	139	57.9
Solyc05g012020.2.1	RIN	MADS-box	MELO3C026300	MADS-L2*	1.29E-99	290	59.1
			MELO3C022316	MADS-L3*	3.04E-88	260	57.0
			MELO3C002049	MADS-L6*	3.02E-79	238	69.5
			MELO3C022516	MADS-L9*	2.76E-67	207	47.0
			MELO3C002050	MADS-L1*	1.10E-60	191	44.1
Solyc03g114840.2.1	MADS1	MADS-box	MELO3C026300	MADS-L2*	1.06E-106	308	63.5
			MELO3C022316	MADS-L3*	4.14E-95	278	59.1
			MELO3C002049	MADS-L6*	8.95E-89	263	56.2
			MELO3C022516	MADS-L9*	1.62E-67	208	50.0
			MELO3C002050	MADS-L1*	1.41E-63	198	46.2
Solyc07g055920.2.1	TAGL1	MADS-box	MELO3C002691	MADS-L4	1.11E-94	278	64.6
			MELO3C007181	MADS-L5	2.07E-94	277	64.3
			MELO3C022209	MADS-L7	2.72E-85	253	59.0
			MELO3C022516	MADS-L9*	2.03E-49	162	41.3
			MELO3C018049	MADS-L10	1.30E-47	156	48.6
Solyc06g069430.2.1	FUL1	MADS-box	MELO3C002050	MADS-L1*	3.27E-109	314	64.9
			MELO3C011409	MADS-L8*	1.12E-74	223	74.5
			MELO3C022316	MADS-L3*	2.82E-60	189	45.7
			MELO3C026300	MADS-L2*	2.66E-56	177	54.1
			MELO3C022516	MADS-L9*	4.45E-56	177	55.7
Solyc03g114830.2.1	FUL2	MADS-box	MELO3C002050	MADS-L1*	1.02E-115	330	68.4
			MELO3C011409	MADS-L8*	5.39E-74	221	71.8
			MELO3C026300	MADS-L2*	9.51E-59	186	43.9
			MELO3C022316	MADS-L3*	1.11E-57	183	54.9
			MELO3C002049	MADS-L6*	2.07E-54	175	52.0
Solyc10g006880.2.1	NOR	NAC	MELO3C016540	NAC-L2*	5.59E-121	352	55.0
			MELO3C018237	NAC-L4	6.41E-100	297	49.9
			MELO3C002573	NAC-L5	6.30E-88	268	44.7
			MELO3C012114	NAC-L6	3.26E-78	241	64.7
			MELO3C018242	NAC-L7*	5.04E-78	241	61.9
Solyc11g017470.1.1	NAC4	NAC	MELO3C010632	NAC-L1	3.76E-128	366	62.5
			MELO3C023195	NAC-L3	8.32E-116	335	59.5
			MELO3C018242	NAC-L7*	4.56E-75	231	65.9
			MELO3C016536	NAC-L8	4.94E-75	231	66.5
			MELO3C016540	NAC-L2*	3.79E-73	228	64.1
Solyc02g077920.2.1	CNR	Squamosa-like	MELO3C025597	CNR-L1	1.37E-44	144	60.2
			MELO3C007273	CNR-L2	1.43E-43	139	56.5
			MELO3C002618	CNR-L3	2.42E-41	134	55.7
			MELO3C009639	CNR-L4	3.82E-36	121	66.7
			MELO3C005966	CNR-L5	7.73E-36	125	55.8

Possible fruit ripening-related melon genes were identified in the melon genome reference CM3.5.1 by BLASTp using tomato fruit ripening-related genes [E8, AP2a, HB1, RIPENING-INHIBITOR (RIN), MADS1, TAGL1, FUL1, FUL2, NON-RIPENING (NOR), NAC4 and CNR] as queries. The top five best hits were identified and listed for each tomato gene.

Some MADS-box and NAC transcription factor genes overlapped between queries (asterisks). These melon genes were subjected to Melonet-DB ‘Coexpression viewer’ as shown in Fig. 6.

Fig. 6

Coexpression network analysis of fruit ripening-related homologs by Melonet-DB’s ‘Coexpression viewer’. (A) Coexpression network of fruit ripening-associated candidate genes (Table 2) and volatile-related genes (Fig. 5) was analyzed by Melonet-DB ‘Coexpression viewer’ (http://melonet-db.agbi.tsukuba.ac.jp/cgi-bin/melcor.cgi). Nodes (genes) are connected by edges when their expression patterns are similar for a defined threshold. (B) Coexpression network of fruit ripening-associated genes identified in (A) (red box with dotted line) and other NAC/MADS-box/homeobox transcription factor genes were subjected to Melonet-DB ‘Coexpression viewer’ to identify additional fruit ripening-related genes. The fruit ripening-related cluster (left, magnified coexpression network) was distinguished from other clusters (right). (A, B) Analysis conditions: weight cut-off = 0.03, RSQ cut-off = 0.3, coexpression database = ‘Harukei3 29 different tissues or fruit stages without callus’.

Fig. 7

Digital expression analysis of fruit ripening-associated genes by Melonet-DB’s ‘Gene expression map viewer’. Digital expression analysis of fruit ripening-related melon genes identified in Fig. 6. (A) Gene expression patterns of fruit ripening-associated melon homologs for tomato E8, NON-RIPENING (NOR), RIPENING INHIBITOR (RIN), MADS1, FUL1/2, NAC4 and HB-1. (B) Gene expression patterns of early fruit development-associated melon homologs for tomato CNR and TAGL1. (C) Expression pattern of melon MADS-box like transcription factor genes which exhibited strong expression at DAF 43 and 50.

Finally, we queried for additional MADS-box transcription factor genes, NAC domain transcription factor genes and homeobox transcription factor genes using the BLASTp search tool, and conducted coexpression analysis again to identify further fruit ripening-associated transcription factors in the melon genome. In total, 153 melon genes were identified as putative MADS-box, NAC domain or homeobox transcription factors from the BLASTp query against the Arabidopsis TAIR10 annotation dataset (Supplementary Dataset S3). Coexpression analysis indicated that the fruit ripening-associated cluster was distinguished from other network clusters (Fig. 6B). In addition, some new transcription factor genes were identified as fruit ripening-associated genes (MELO3C022921, MELO3C018088, MELO3C011979, MELO3C014519, MELO3C013727, MELO3C010678, MELO3C009873, MELO3C002628 and MELO3C015239). Interestingly, two homeobox transcription factor genes (MELO3C022921 and MELO3C018088) show strong association with CmAAT1 and CmAAT2 (Fig. 6B). As the expression of CmAAT1 and CmAAT2 was strongly induced during the late stage of fruit development (DAF 43 and 50; Fig. 5), but was absent around DAF 22–36, they are located at marginal positions in the fruit ripening-associated network (Fig. 6B). In agreement with the coexpression status, the expression of MELO3C022921 and MELO3C018088 was found to be induced at DAF 43 and 50 (Fig. 7C), suggesting their role in the late stage of fruit ripening.

In conclusion, we demonstrated that our tools “‘Gene expression map viewer”’ and “‘Coexpression viewer”’ are effective to identify computationally fruit ripening-associated melon genes. It is expected that the in silico analysis strategy shown in this study can be applied to the study of other genes or biological phenomena. However, the current gene expression dataset in Melonet-DB does not contain some important tissues or other fruit ripening stages (i.e. germinating seed and post-harvest fruit). Given that ethylene emission was observed during post-harvest stages (Supplementary Fig. S1), it appears valuable to include post-harvest fruits in the coexpression dataset to explore the relationship between fruit ripening-associated transcription factors (Fig. 6) and ethylene biosynthesis. In addition, it is also important to generate ethylene biosynthesis mutants and compare global gene expression profiles between the wild type and such mutants. This will contribute to our understanding of gene expression networks in the light of hormonal regulation (i.e. ethylene-dependent and independent regulation). In addition to the Melonet-DB database, we are also developing an artificially induced mutant population as well as genome editing methodology in the ‘Harukei-3 melon. By combining such genetic experiments, it will be possible to clarify whether genes identified in this study are actually involved in fruit ripening or other process(es). Such efforts will contribute to the functional genomics study of melon and fleshy fruit plant species.

Materials and Methods

Plant materials and growth conditions

Seeds of the ‘Earl’s favorite Harukei-3’ melon were obtained from the National Agriculture Research Organization (NARO) Genebank. Seeds were germinated in soil composed of ‘Na-tera’ (Mitsubishi Chemical Agri Dream) in the dark at 25–28 °C in April or May, and seedlings were grown under 9 h light (25 °C)/15 h dark (20 °C) conditions until true leaf expansion. Plants were transferred to a soil composed of ‘Coco-bag’ (Toyotane) and grown under greenhouse conditions. Plants were irrigated with a hydroponic solution composed of ‘Tank-mix A&B’ (OAT Agrio). Female flowers were hand pollinated with the pollen from male anthers that were obtained from the same plant (self-pollination), and one or two fruits were retained in each plant. In the case of disease emergence (i.e. powdery mildew, whitefly or canker), chemicals were applied to the plant. For the isolation of total RNA, harvested plant tissues were frozen in liquid nitrogen and stored at − 80 °C until use.

Isolation of total RNA and RNA-seq data acquisition

For isolation of total RNA, frozen plant tissues were ground with a multi-bead shocker (Yasui Kikai) that had been pre-cooled with liquid nitrogen. Total RNA was extracted and purified using a Maxwell® 16 LEV Plant RNA Kit (Promega) according to the manufacturer’s protocol. Isolated total RNA was stored at − 80 °C until use. RNA-seq was performed using an Illumina Hiseq® 2500 with 125 bp paired-end mode (Eurofins Genomics). The DRA accession number for these RNA-seq data is DRA006228.

RNA-seq data analysis

RNA-seq data were first subjected to quality filtering by the fastx-toolkit with parameters “‘-Q 20 -P90”’, and then NGS adaptor sequences were removed if they were present at the terminal end of short reads. To eliminate possible contaminated fungal/bacterial ribosomal reads that might originate from microorganisms attached to the melon plants grown in the greenhouse, short reads were aligned to the Greengene 16S ribosomal sequence dataset (McDonald et al. 2012) and the fungal rDNA internal transcribed spacer (ITS) sequences dataset (Schoch et al. 2012) using bowtie2 software (Langmead and Salzberg 2012). Unmapped reads were extracted and used for subsequent analysis of melon gene (total 34 Gb, average 4.6 million reads per sample). The average rate of reads that were aligned to the bacterial/fungal datasets was 0.2% (total 561,777 reads). Gene expression levels were calculated as FPKM values by the Tophat-Cufflinks pipeline with “‘very-sensitive”’ mode (Trapnell et al. 2010, Kim et al. 2013) using the CM3.5.1 melon nuclear genome sequence reference that was combined with chloroplast and mitochondrial genome sequence information (Rodríguez-Moreno et al. 2011, Garcia-Mas et al. 2012). Of the 27,427 genes present in the CM3.5.1 melon genome reference, 20,752 genes had FPKM with >1.0 in at least one melon tissue, and the remaining 6,675 genes were removed as non-expressing genes. FPKM values were transformed into relative values (maximum value as 1.0 in each gene) after normalization by two melon genes, MELO3C026304 and MELO3C011153, which are homologs of the constitutively expressed Arabidopsis CAP-binding protein gene (AT5G44200.1) or ubiquitin-conjugating enzyme gene (AT5G25760.1), respectively. The gene expression dataset was subjected to WGCNA (ver. 1.5.1) implemented in R software (ver. 3.2.3) (Langfelder and Horvath 2008, Langfelder and Horvath 2012) to conduct co-expression analysis. In PCA, log₂-transformed values were subjected to the ‘princomp’ function implemented in R software. Graph images showing the gene expression pattern in 30 melon tissues (Fig. 3B) were generated by perl scripts.

InterProScan search and ID enrichment analysis

InterProScan search was conducted using the deduced protein sequence information of 34,848 transcripts in the CM3.5.1 melon reference (Jones et al. 2014). Of these, 23,340 transcripts were found to encode protein with an IPR ID (Supplementary Dataset S1). ID enrichment was analyzed in co-expression groups identified by WGCNA by perl scripts (Supplementary Dataset S2).

Identification of putative MADS-box, NAC and homeobox transcription factor genes

To identify melon homologs for tomato fruit ripening regulators (E8, RIN, MADS1, TAGL1, FUL1/2, NOR, NAC4 and CNR), a BLASTp search was conducted using the deduced protein sequences of tomato genes as queries and the entire protein sequence dataset of the CM3.5.1 melon genome reference as database. Because E-values varied significantly and were dependent on tomato queries, the top five melon gene hits were identified and subjected to co-expression analysis implemented in Melonet-DB. To identify putative MADS-box, NAC and homeobox transcription factor genes in the melon genome, a BLASTp search was conducted using the entire protein sequence dataset in the CM3.5.1 genome reference as queries and Arabidopsis TAIR10 protein sequence as the database. Based on gene annotation information of Arabidopsis genes, putative MADS-box, NAC and homeobox transcription factor genes were identified in melon by setting thresholds as follows: E-value = 1e^–50, % similarity = 50 (Supplementary Dataset S3).

Development of web-application tools in the Melonet-DB

“‘Gene expression map viewer”’ and “‘Coexpression viewer”’ tools were programmed mostly in perl and installed in a LAMP server located within the University of Tsukuba (http://melonet-db.agbi.tsukuba.ac.jp/). The latter tool utilizes the Cytoscape web for visualization of gene to gene correlation networks (Lopes et al. 2010). Although databases are currently available for 30 melon tissues in the Melonet-DB website, new datasets will be available when we conduct additional RNA-seq studies. Basic usage instructions of our tools are also available on the Melonet-DB website (http://melonet-db.agbi.tsukuba.ac.jp/cgi-bin/help.cgi).

Supplementary Data

Supplementary data are available at PCP online.

Funding

This work was supported by the program of the SIP project [(Cross-ministerial Strategic Innovation Promotion Program, the Cabinet Office) to S.N and H.E.]; the Japan Society for the Promotion of Science [Grant-in-Aid for Young Scientists B to R.Y.]; and the Japan Science and Technology Agency (JST) [PRESTO programme to R.Y.].

Acknowledgments

We thank Drs. Toru Ariizumi and Naoya Fukuda (University of Tsukuba) for valuable and helpful discussions, and Ms. Chiho Mito for her technical assistance.

Disclosures

The authors have no conflicts of interest to declare.

Abbreviations

Abbreviations
 
• AAT

  alcohol acyl-transferase

 
• ADH

  alcohol dehydrogenase

 
• CG

  coexpression group

 
• CNR

  COLORLESS NON-RIPENING

 
• DAF

  day after fertilization

 
• FPKM

  fragments per kilobase of transcript per million fragments sequenced

 
• NGS

  next-generation sequencing

 
• NOR

  NON-RIPENING

 
• PCA

  principal component analysis

 
• QTL

  quqntitative trait locus

 
• RIN

  RIPENING INHIBITOR

 
• RNA-seq

  RNA-sequencing

 
• WGCNA

  weighted gene correlation network analysis