Transcriptome analysis of the oriental melon (Cucumis melo L. var. makuwa) during fruit development

Plant materials

Oriental melon plants (Gotgam) were grown in a greenhouse with plastic pots containing a soil mixture. The temperature was maintained at 28 °C during the day and 24 °C during the night Supplemental lights were used to provide an 18-h light period. Flowers were hand pollinated and at least three fruits were harvested from different plants at 8, 16, 24 and 32 DAP. Tissues were immediately frozen in liquid nitrogen and stored at −80 °C until analysis.

Total RNA extraction, de novo assembly and annotation

Total RNA was isolated from each sample of fruit flesh using TRIzol reagent according to the manufacturer’s instructions. Purified total RNAs were used to synthesize cDNAs, which were amplified according to the Illumina RNA-Seq protocol and sequenced using an Illumina HiSeq2000 system, producing 17.9 Gb 101-bp paired-end reads in total from one sequencing run with pooled by three independent fruits in each stage. The sequence data were filtered by quality score (Q ≥ 20) and minimum length (≥25 bp) using SolexaQA (Cox, Peterson & Biggs, 2010). The quality-checked sequence reads from different tissue samples were assembled de novo using two software tools based on de Bruijn graphs: Velvet (v1.2.07) (Zerbino & Birney, 2008) and Oases (v0.2.08) (Schulz et al., 2012). The assembled transcripts were used for DEG extraction or functional annotation by edgeR or BLASTP with Cucumis melo or Arabidopsis thaliana as queries (Kim et al., 2015; Klukas & Schreiber, 2007). The annotated genes were used for further analyses including gene ontology (GO) and KEGG analyses.

Differentially expressed genes and clustering analysis

DEGs were selected by comparing two distinct stages of fruit development. The reads with each sequence tag were mapped to the assembled loci using Bowtie (mismatch ≤2 bp); the number of cleanly mapped reads for each locus was determined; and the data were normalized using the DESeq library in R (Anders & Huber, 2010). The final DEGs were identified by a change in read coverage of at least twofold, a false discovery rate (FDR) ≤ 0.05, and a read count ≥ 100 in both samples. The clustering of the DEGs was performed using the Mfuzz program (Futschik & Carlisle, 2005; Kumar & Futschik, 2007).

Functional enrichment and gene ontology analyses

The assembled loci were annotated to the GO database (downloaded from http://www.geneontology.org/) using BLASTP (e-value = 1e⁻⁰⁶). The annotation was carried out using the GO classification results from the Map2Slim.pl script (Kim et al., 2016a). Functional enrichment analysis was carried out using DAVID, a web-accessible program that provides a comprehensive set of functional annotation tools to infer biological meaning from large lists of genes (Huang Da, Sherman & Lempicki, 2009; Huang, Sherman & Lempicki, 2008). Fisher’s exact test was used to analyze the gene lists annotated with the TAIR identifications of the transcripts with respect to the GO terms under the following criteria: count ≥10; FDR ≤0.05.

Kyoto Encyclopedia of Genes and Genomes analysis

KEGG analysis was performed to investigate the biochemical pathways involved in fruit development using the DEGs among the KM transcripts. A. thaliana orthologs of the KM DEGs in each cluster were identified, and KEGG analysis was performed with the Plant GSEA platform (http://structuralbiology.cau.edu.cn/PlantGSEA/analysis.php) (Yi, Du & Su, 2013) using the default parameters.

De novo sequence assembly and functional annotation

The KM fruit has a rounded, rectangular shape, medium size, deep-green skin, and light-red flesh. The progression of growth and development in the KM fruit is highly reproducible and includes visible external and internal morphological changes. The fruit was still small at 8 DAP. From 8 DAP to 16 DAP, the diameter and length of the fruit increased rapidly. After 16 DAP, the fruit volume continuously increased to the maximum, and the color changed to dark green (Fig. 1A).

For transcriptome analysis of the fruit development and maturation, we collected samples at two developmental stages (8 DAP and 16 DAP) and two maturation stages (24 DAP and 32 DAP) for total RNA extraction. Sequencing with Illumina HiSeq 2000 produced 17.9 Gb of paired-end data in total. To determine a unigene set, we performed de novo transcriptome assembly with Velvet/Oasis and assembled 100,875 transcripts. In general, the assembled transcripts contain the alternative spliced form of transcripts or partial assembled transcripts. To remove these transcripts and determine representative loci, we performed reciprocal blast and the longest and non-redundant transcripts were determined. Then, BLASTx was performed with remained transcripts using Non-Redundant (NR) database and 1,722 of transposable element related genes were removed. Finally, 47,666 transcripts were remained as representative loci (Table 1 and Fig. S1). We then annotated the representative loci by BLASTP (e-value ≤ 1e⁻⁰⁸) against melon (CM V3.5) and Arabidopsis genes (TAIR 10), as reported previously (Kim et al., 2016a). Furthermore, we also performed ISGAP pipeline to determine more precise transcripts by removing intron contaminants (Kim et al., 2015). Of the representative loci, 35,364 were assigned to melon reference proteins, and 33,963 were assigned to Arabidopsis proteins. For further analyses, we used the representative loci, including the functionally annotated loci, to remove alternative transcripts and assembly errors.

Gene ontology classification of differentially expressed genes

To examine differential gene expression, we mapped the trimmed raw reads to the representative loci to investigate KM specific transcripts and their expression patterns and determined the genes that were expressed differentially between the indicated time points (Fig. 1B). The results showed that 253 and 446 genes were up-regulated and down-regulated, respectively, at 16 DAP relative to 8 DAP. In addition, at 24 DAP and 32 DAP, 792 and 2,794 genes were up-regulated, respectively, while 1,523 and 3,034 genes were down-regulated, respectively. These results showed that many genes were transcriptionally regulated during fruit maturation and that various metabolic pathways involved in fruit development or maturation were dynamically regulated at the transcriptional level. To investigate the functions of the DEGs in each developmental stage, we performed a GO analysis to classify the DEGs (Fig. 2 and Fig. S2). We assigned a total of 4,557 DEGs to 38 functional categories: 19 ‘biological processes,’ 11 ‘cellular components,’ and eight ‘molecular functions’ (Fig. 2). In each of those categories, ‘cellular process’ (7.75%), ‘cell and cell part’ (20.84%), and ‘binding and catalytic activity’ (10.60%) were the most abundant GO terms, respectively. The GO analysis revealed that both the up-regulated gene set and the down-regulated gene set were significantly enriched with genes involved in ‘catalytic activity’ and ‘binding activity’. Those results suggest that genes involved in catalytic pathways were regulated dynamically during fruit maturation.

Clustering and Kyoto Encyclopedia of Genes and Genomes pathway enrichment of the differentially expressed genes

To better understand the transcriptional regulation of the gene expression involved in fruit development and maturation, we performed a cluster analysis of the DEGs (Futschik & Carlisle, 2005; Kumar & Futschik, 2007). We classified a total of 5,173 DEGs into 14 types of clusters based on the modulation of expression patterns (Fig. 3 and Table 2).

To investigate the biological functions of the candidate genes potentially involved in fruit development and maturation, we used the 5,173 DEGs for further KEGG analysis (Kanehisa et al., 2012). We identified orthologs of the KM loci in A. thaliana and applied the orthologs for KEGG analysis. The KEGG analysis identified eight pathways including starch and sucrose metabolism (ath00500), carotenoid biosynthesis (ath00906), phenylpropanoid biosynthesis (ath00940), and pentose and glucoronate interconversions (ath00040) (Table 3). Those results indicated that genes involved in metabolic pathways were transcriptionally regulated during fruit development (8 DAP and 16 DAP) and maturation (24 DAP and 32 DAP). We mapped a total of 48 DEGs to starch and sucrose metabolism pathways and grouped them into 11 clusters (Fig. 4). Among various pathways, glucose metabolic pathways showed finely tuned expression patterns of genes. Two glucose metabolic pathways showed high gene expression patterns in both fruit development and maturation stages (red circles in Fig. 4A) and one pathways showed high gene expression patterns in maturation stage (orange circle in Fig. 4A). Six loci (Locus_27541, Locus_1366, Locus_1187, Locus_7852, Locus_30074, and Locus_8322) were involved in three glucose-metabolism pathways with different substrates such as UDP-glucose, cellulose, and maltose (Fig. 4A). These loci were highly expressed during fruit development and subsequently down-regulated during fruit maturation. In contrast, three loci (Locus_7857, Locus 12354, and Locus_5309) showed increased expression during fruit maturation. Those results indicated that nine genes in three glucose-metabolism pathways were involved in the regulation of glucose contents in the oriental melon fruit. Genes involved in fructose metabolism were highly expressed during fruit development (8 DAP) and subsequently down-regulated during fruit maturation (purple circle in Fig. 4A). Furthermore, genes involved in UDP-glucose and sucrose-6P metabolism, both precursors of sucrose, were up-regulated at 8 DAP and down-regulated during fruit maturation (purple circles in Fig. 4A).

The KEGG analysis identified 13 genes involved in carotenoid biosynthesis. In the carotenoid-biosynthesis pathway, phytoene synthase (PSY, EC 2.5.1.32) (Bouvier, Rahier & Camara, 2005), an enzyme that catalyzes the synthesis of 40-carbon phytoene from the condensation of two geranylgeranyl diphosphate (GGPP) molecules (Bouvier et al., 2005), was up-regulated during fruit maturation. This step is the first step in carotenoid biosynthesis. In contrast to PSY, three carotenoid hydroxylases (CYP97A3) (Kim & DellaPenna, 2006) were up-regulated during fruit development and subsequently down-regulated during fruit maturation. For xanthophyll biosynthesis, three 9-cis-epoxycarotenoid dioxygenases (Schwartz et al., 1997) (NCEDs), which yield a C15 intermediate xanthoxin, were differentially expressed during fruit development and maturation. NCED5 was up-regulated during fruit development, whereas NCED3 and NCDE4 were up-regulated during fruit maturation. In addition, three CYP707 genes (Millar et al., 2006) involved in ABA oxidation were up-regulated during fruit development (Fig. 5). In mature melon fruits, beta-carotene is the most abundant carotenoid, with trace amounts detected even in the white-fleshed Piel de sapo (Saladie et al., 2015). Our results corresponded to previous results and indicated that the expression patterns of carotenoid biosynthesis genes were regulated to increase the carotenoid contents during fruit maturation.