Melon (Cucumis melo L.) is an important cultivated cucurbit and produces an aromatic, sweet fruit. Melons exhibit significant variation in terms of botanical phenotypes and fruit types, highlighting the genetic diversity of the Cucurbitaceae family (Mliki et al., 2001). Six melon subspecies are cultivated: cantalupensis, reticulatus, inodorus, acidulous, saccharinus, and makuwa (Liu et al., 2004).

Melon is a diploid species, with a basic number of chromosomes (x = 12 [2x = 2n = 24]) and an estimated genome size of 450 Mb, similar to that of rice (419 Mb). The melon genome is being sequenced as part of the Spanish Genome Initiative (MELOGENOMICS). Moreover, BAC libraries, high-resolution genetic maps, oligo-based microarrays, and a large number of transcriptome sequences (RNA-Seq and expressed sequence tag(EST)) for melon are also available as genetic and genomic tools.

Oriental melon (C. melo L. var. makuwa) is cultivated in large areas of Asia, especially throughout the temperate regions of China, Japan, and Korea. KM (Korean landrace; Gotgam) is one of the major landrace of oriental melon in Korea. This land-race contains more nutrients and has useful traits with greater disease resistance than other varieties. NW is high quality breeding line of NongWoo Bio Company with deep yellow, middle sized, oval type fruit and contained high sugar contents. Although oriental melon is economically valuable in Asia and exhibits genetic features distinct from the other subspecies, few reports of genome-scale studies of oriental melon have been published. Therefore, development of genomic resources and markers of these oriental melon will be important for the characterization of genetic markers linked to highly desirable traits and oriental melon breeding.

Due to the high-throughput capacity of next-generation sequencing (NGS) technology, which was developed in the last decade, transcriptome analysis has become widely used for genome-scale studies. Transcriptome analysis can be used to profile gene expression and identify novel transcripts, splicing isoforms, and sequence variations, including single-nucleotide polymorphisms (SNPs) and simple sequence repeats (SSRs). In the present work, we generated a total of 67,440,566,178 raw sequence reads (67.4 Gb) from the female and male flowers, leaves, roots, and fruit of two oriental melon varieties (KM and NW). A total of 64,998 transcripts from KM and 100,234 from NW specimens were de novo assembled, with transcript N50 values of 939 and 1,138, respectively. A total of 259 SNP and SSR markers were developed, and genotype- and species-specific genes were identified. The transcripts generated in the present study will be a useful resource for the characterization of gene expression patterns and traits in oriental melon.

Plant materials and mRNA sequencing

To characterize the oriental melon transcriptome and increase the sequence coverage of C. melo, sequence libraries for two oriental melon cultivars were generated. Specimens of the KM (Korean landrace) and NW (Breeding line of NongWoo Bio Co.) cultivars were provided by NongWoo Bio Company for use in the study (Supplementary Fig. S1). Thirty two seeds of each cultivar were sowed in tray and grown in a greenhouse until 3?4 true leaf stage. Ten healthy plants of each cultivar were transplanted into plastic pot and grown until flowering stage in greenhouse in which temperature was maintained 28°C in daytime and 24°C at night. Total RNA was isolated from the leaves, roots, female and male flowers, and fruit of both cultivars using TRIzol reagent according to the manufacturer’s instructions. RNA from three biological replicates were pooled before cDNA synthesis. Purified total RNA was used to synthesize cDNAs, which were amplified according to the Illumina RNA-Seq protocol and sequenced using an Illumina HiSeq2000 system, producing 2.6 G 101-bp paired-end reads.

De novo assembly and annotation

Sequence data with a quality score above 20 (Q ≥ 20) were extracted using SolexaQA (Cox et al., 2010; Kim et al., 2014). Sequence reads from different tissue samples were de novo assembled using two software tools based on the de Bruijn graph algorithm: Velvet (v1.2.07)(Zerbino and Birney, 2008) and Oases (v0.2.08)(Schulz et al., 2012). A K-mer value of 59 was considered indicative of the optimum length for de novo assembly of oriental melon sequence reads. A schematic illustration of the process is shown in Fig. 1.

The quality-checked reads from each tissue were merged and used for transcript assembly. The assembled transcripts were validated by direct comparison with gene sequences in the SEEDERS plant annotation database using BLASTX (evalue ≤ 1e^?05) (Altschul et al., 1990). Protein sequences with the highest similarity were retrieved for further analysis. Short reads of KM transcripts were mapped to the MELONOMICS melon genome (http://melonomics.net/) (Garcia-Mas et al., 2012), and CDS regions were defined as in the melon genome.

Functional enrichment analysis

For gene ontology (GO) term analysis, the assembled loci were annotated to the GO database (downloaded from http://www.geneontology.org/) using BLASTP (e-value = 1e^?06). GO term annotation was carried out using GO classification results from the Map2Slim.pl script. Functional enrichment analysis was carried out using DAVID, a web-accessible program that provides a comprehensive set of function annotation tools for inferring biological meaning from large lists of genes (Huang da et al., 2009a; 2009b). Fisher’s exact test was used to analyze the gene lists annotated with the TAIR identifications of the transcripts with respect to GO terms, under the following criteria: counts ≥ 10; false discovery rate (FDR) ≤ 0.05.

Short-read counting and tissue-specificity scoring

Illumina sequencing was used to generate mRNA libraries for the various oriental melon tissues examined. Reads for each sequence tag were mapped to the assembled loci using Bowtie (mismatch ≤ 2 bp); the number of clean mapped reads for each locus was determined, and the data were normalized using the DESeq library in R. Only transcripts with a tag count ≥ 50 were retained for further analysis. Genes differentially expressed between samples were identified based on the fold-change in expression, with the results analyzed by t-test. The FDR determined from multiple tests and analyses was applied to calculate the p-value threshold via DESeq. For reliable determination of tissue specificity, Fisher’s exact test was used to compare the proportion of a given transcript among all transcripts in each tissue examined. Tissue-specific genes were selected based on read counts of > 50 in the target tissue and < 10 in the other tissues. The degree of homology between the sequences of KM, NW, and melon transcripts was determined using BLASTn/tBLASTn (e-value ≤ 1e^?100; identity ≥ 90).

Identification of SNPs

To identify SNPs, a quality check of the KM and NW raw reads was performed using the Solexa QA package. The raw reads were aligned against melon mRNA sequences using TopHat, with modified default parameters (mismatches [?N] = 1; maximum insertion length = 1; minimum intron length [?i] = 50; maximum intron length [?I] = 14,018; mate inner distance [?r] = 350; segment mismatches = 1; maximum segment intron = 100), and the results were saved as a BAM file for further analysis using SAMTools (Kim et al., 2014; Li et al., 2009).

Using the varFilter command in SAMTools, SNPs were called only for variable positions with a minimum mapping quality (?Q) of 30. The minimum and maximum read depths were set at 3 and 1,000, respectively. Significant SNP sites among the sequences of transcripts from KM, NW, and melon were identified using a Perl script developed in-house (Supplementary Fig. S2).

Identification of SSRs

To identify SSRs, assembled transcripts of NW specimens were formatted according to the SSR Locator’s protocol (da Maia et al., 2008). Perfect SSRs (designated ‘P-type’ SSRs) forming dimer to hexamer motifs with more than five repeat units and located more than 100 bp from other SSRs were selected. Imperfect SSRs (designated ‘I-type’ SSRs) were selected by allowing for 5 bp of erroneous sequence. Previously reported criteria (Garg et al., 2011; Kong et al., 2007) were used to select SSRs. Primers were designed using Primer 3 in the SSR Locator. Using the designed primer sets, virtual PCR was performed with SSRs from NW specimens according to the SSR Locator’s protocol. Transcripts containing one P-type SSR were selected and used for the development of KM, NW, and melon markers. The following selection criteria were used for the primers: (i) the expected amplicon size should be the same as that of the virtual PCR; and (ii) the primer sets for the different SSRs should not overlap. The selected primer sets were used for virtual PCR analysis of the KM and melon sequence data to distinguish NW-specific marker candidates (Supplementary Fig. S3).

Genetic linkage map construction

An oriental melon genetic map was constructed using MAPMAKER 3.0/EXP (Lander et al., 1987) with 234 dCAPS (Neff et al., 1998) and 25 SSR markers. F2 population NW and KM specimens were used for mapping. Recombination fractions were converted to map distances in centimorgans (cM) using the Kosambi mapping function (Kosambi, 1943).

De novo assembly of the oriental melon transcriptome

To perform transcriptome analysis, RNA-Seq data were generated from five different tissues (female and male flowers, fruit, leaves, and roots) of two oriental melon cultivars (KM and NW). In total, 30.5 Gb (251,752,490 raw reads) and 36.8 Gb (287,233,170 raw reads) of KM and NW sequence data, respectively, were generated using an Illumina HighSeq 2000 (Supplementary Table S1). The quality of the sequence data (Q ≥ 20) was assessed using SolexaQA, and the reads were trimmed and sorted by length using the DynamicTrim and LengthSort programs, respectively.

Transcripts of each oriental melon cultivar were assembled using Velvet (v1.2.07) and Oases (v0.2.08) (k-mer = 59), based on de Bruijn graphs. These transcripts were used to construct extended transcripts using Velvet, followed by Oases (Fig. 1). A total of 64,998 KM transcripts were generated, with a mean length of 706 bp and N50 length of 939 bp; the transcript length ranged from 200 to 13,444 bp. For the NW cultivar, 100,234 transcripts were assembled, with a mean length of 739 bp and N50 length of 1,138 bp; the transcript length ranged from 117 to 11,659 bp. The transcripts from each cultivar were then clustered, resulting in 49,409 predicted loci for KM and 51,557 loci for NW (Table 1).

Functional annotation and classification of oriental melon transcripts

Putative functions of the assembled transcripts were annotated using BLASTP (e-value ≤ 1e^?06) with the SEEDERS non-redundant protein database. Of the 64,998 KM transcripts, 36,871 were assigned to 21,363 reference proteins, and 64,149 of 100,234 NW transcripts were assigned to 21,914 reference proteins (Table 2). To classify the functions of the assembled oriental melon loci, GO term analysis was performed using TAIR identification information. A total of 42,386 KM and 42,743 NW transcripts were assigned to 23 functional categories: 13 ‘biological process’ categories, 7 ‘cellular component” categories, and 3 ‘molecular function’ categories (Fig. 2). For both the KM and NW transcripts, ‘cellular process’, ‘cell and cell part’, and ‘catalytic activity’ were the most common terms in the ‘biological process,’ ‘cellular component,’ and ‘molecular function’ categories, respectively. These GO term data will be used for further studies of the characteristics of oriental melon by functional profiling, prediction of gene function, and functional categorization of genes (Rhee et al., 2008).

Expression of tissue-specific locus candidates

Constitutive promoters, such as those for ubiquitin and 35S, are used in plant genetic engineering to express genes of interest in a wide range of species (Brisson, 1984; Cornejo et al., 1993). However, overexpression using constitutive promoters may lead to undesirable pleiotropic effects in transgenic plants (Hsieh et al., 2002; Kasuga et al., 1999). The use of tissue-specific promoters with particular developmental expression patterns has been suggested as a strategy to avoid such undesirable pleiotropic effects (Kasuga et al., 2004). Therefore, the development of tissue-specific promoters capable of driving transgene expression is an important area of research in plant genetic engineering (Potenza, 2004).

For unbiased detection of tissue-specific expressed transcripts, statistical analysis of a large number of raw reads was performed. Fisher’s exact test was used to compare the proportion of given transcripts among all transcripts in the different tissues. Non-normalized oriental melon cDNA libraries were prepared from the female and male flowers, fruit, leaves, and roots of each of the oriental melon cultivars. The sequence tags of short reads were mapped to each transcript using Bowtie, and the number of mapped reads for each transcript was then determined. The mapped reads were normalized using the DESeq library in R script. Statistical analyses identified 1,169 and 2,504 tissue-specific transcripts in the KM and NW cultivars, respectively. For the KM cultivar, 144, 107, 231, 256, and 431 transcripts were specific to the female flower, male flower, fruit, leaves, and roots, respectively. In the case of the NW cultivar, 75, 1,121, 368, 150, and 790 transcripts were specific to the female flower, male flower, fruit, leaves, and roots, respectively (Supplementary Tables S2 and S5). The functions of the KM and NW transcripts were predicted by identifying orthologues using melon mRNA sequences. In total, 16,873 of 27,427 melon mRNAs were identified as orthologues of 48,144 KM transcripts and 66,216 NW transcripts. Among these transcripts, 914 were KM-specific, identified based on 746 melon mRNAs, whereas 1,070 NW-specific transcripts were identified based on 573 melon mRNAs (Supplementary Table S3). These tissue-specific candidates will be validated using RT-PCR, and the promoter regions will be investigated for cloning tissue-specific promoters. Furthermore, functional studies of tissue-specific genes will provide additional insights into plant development.

Identification of SNPs and SSRs

Molecular markers are important resources for constructing high-density genetic maps such as those used in crop breeding and for the identification of traits of interest. Since NGS technology was developed, many plant genomes have been sequenced, including that of melon (Garcia-Mas et al., 2012). In addition, a large amount of sequence data for melon has been accumulated over the past several years (Gonzalez-Ibeas et al., 2007; 2010; Portnoy et al., 2011; Rodriguez-Moreno et al., 2011). SSRs and SNPs are increasingly used in the construction of melon genetic maps (Blanca et al., 2011; 2012; Diaz et al., 2011; Kong et al., 2011). SNPs and SSRs were identified among the KM, NW, and melon transcripts using the assembled transcripts and melon mRNA sequences. A total of 7,871 SNPs covering 2,156 loci and 3,110 transcripts were identified between the KM and NW cultivars (KM/NW), and 4,752 SNPs were identified in exon regions. Between the KM and melon sequences (KM/melon), 3,730 SNPs were identified covering 1,063 loci and 1,547 transcripts, and 2,297 SNPs were identified in exons (Table 3; Supplementary Table S4).

The distribution of synonymous and non-synonymous SNPs among the 12 melon chromosomes was also investigated. The number of synonymous and non-synonymous SNPs in chromosome 2 was significantly larger in KM/NW and KM/melon (Fig. 3). The frequency of SNP occurrence between KM and NW was expected to be low, as the sequences were derived from two near-isogenic lines. However, the number of SNPs between KM/NW was larger than between KM/melon. The melon samples used for genome sequencing were double-haploid line, derived from the cross between PI 161375 (Song-whan Charmi) (SC) and the ‘Piel de sapo’ (PS) (Oliver et al., 2001). The NW line was bred from the cross between EunCheon type commercial Fl variety and Chinese landrace melon. Furthermore, EunCheon is derived from the cross between Japanese landrace Charmi and small melon. Half of genome sequences in melon reference were Songwhan Charmi (PI 161375) while NW genome sequence were consisted with lots of melon landrace genome. Thus genome variations of KM and NW were higher than those of KM and melon reference. Consequently, fewer SNPs were identified in KM/melon compared with KM/NW. Data regarding the synonymous and non-synonymous SNPs of NW and melon compared with KM are provided in Supplementary Tables S6 and S7. The GO terms for the KM and NW transcripts or KM and melon mRNA sequences with SNPs were sorted, and the rate of synonymous and non-synonymous SNPs was determined (Supplementary Table S8).

dCAPS primers were designed based on 277 SNPs of KM and NW transcripts and used to screen polymorphic markers. dCAPS primer sets for 245 of the 277 SNPs were moderately amplified in both oriental melon cultivars. The amplified products were digested using restriction enzymes specific to sites in the primer sequences. A total of 234 PCR products exhibited polymorphism between the KM and NW cultivars, and 16 SNP markers are shown in Supplementary Fig. S4.

A general screening of the KM and NW transcript dataset and melon mRNA sequences for unigenes was performed for the presence of di-, tri-, tetra-, penta-, and hexa-SSR motifs. Motifs ranging from dimers to hexamers with more than five repeat units were selected (Garg et al., 2011; Kong et al., 2007). A total of 10,709 SSRs were identified among 9,938 KM transcripts (Supplementary Table S9). The major types of SSRs identified were dinucleotides (6,429), followed by trinucleotides (3,736), tetranucleotides (318), pentanucleotides (127), and hexanucleotides (99) (Table 4). The most frequent SSR motif was GA/TC (1,776), followed by AG/CT (1,574), AT/AT (1,171), TA/TA (1,148), and GAA/TTC (778) (Supplementary Table S10). In NW, 15,662 SSRs were identified among 14,436 transcripts (Supplementary Table S9). The major SSR motifs were dinucleotides (8,708) and trinucleotides (6,219), followed by tetra-nucleotides (412), hexanucleotides (166), and pentanucleotides (157) (Table 4). The GA/TC motif exhibited the highest frequency (2,731), followed by AG/CT (2,401), GAA/TTC (1,399), and AT/AT (1,310) (Supplementary Table S10).

For NW, primers were designed for SSR marker candidates from the transcripts. Transcripts containing one SSR marker candidate were selected for screening polymorphisms among KM, NW, and melon using virtual PCR (Supplementary Table S11). The presence and size of amplicons resulting from the virtual PCR analyses were compared between NW and KM or NW and melon (Supplementary Table S12) for use in further marker development. Of 8,052 SSR marker candidates, 64 were selected for PCR analysis to determine the presence of polymorphisms between KM and NW. A total of 25 SSR markers exhibited a polymorphic pattern between the two cultivars; the results for 16 markers are shown in Supplementary Fig. S5.

The polymorphisms of the 64 SSR markers and 277 dCAPS markers designed from SNP markers were analyzed in the parent lines, and 25 SSR and 234 SNP markers exhibited polymorphisms. All of these polymorphic markers were screened in 94 F2 population plants and exhibited co-dominant type Mendelian segregation. The genotypes of these 259 markers were used to construct an oriental melon genetic linkage map consisting of 12 linkage groups covering 926 cM, with an average map distance between markers of 3.7 cM. A total of 248 markers were assigned to the 12 linkage groups, and 11 markers were unlinked (Fig. 4). The chromosome numbers for all of the linkage groups were determined based on alignment of our DNA marker sequences with the melon genome sequence data (Garcia-Mas et al., 2012) regarding the markers are provided in Supplementary Table S13.

The oriental melon genetic linkage constructed based on the 25 SSR and 234 dCAPS SNP markers spanned 926 cM. Although not all of the markers were distributed evenly among the 12 linkage groups, analysis using the Χ² goodness-of-fit test showed no significant distortions from the expected Mendelian ratio for any of the markers. The largest linkage group was on Ch6, spanning 114 cM, whereas the smallest was on Ch11, spanning 28.2 cM. Because no markers common to other previously published melon linkage maps were used in our study, it was difficult to estimate the consensus between the genome structures. A previously reported integrated melon genetic map was constructed using 1,592 markers from 8 independent mapping experiments and spanned 1,150 cM across the 12 linkage groups (Diaz et al., 2011). In this integrated map, the genetic length of the linkage groups ranged from 73 to 119 cM. The significantly shorter genetic length on Ch11 and lower resolution of linkage groups (8 gaps of more than 20 cM) in our map compared with the previously reported integrated map are most likely due to our development of marker types only from transcriptome data or to the use of an insufficient number of markers. It will therefore be necessary to use more markers derived from the genome sequence database and consensus markers positioned based on several mapping experiments to construct a linkage map in which the markers are evenly distributed and cover all of the genome.