Comparative genomics of muskmelon reveals a potential role for retrotransposons in the modification of gene expression

Melon exhibits substantial natural variation especially in fruit ripening physiology, including both climacteric (ethylene-producing) and non-climacteric types. However, genomic mechanisms underlying such variation are not yet fully understood. Here, we report an Oxford Nanopore-based high-grade genome reference in the semi-climacteric cultivar Harukei-3 (378 Mb + 33,829 protein-coding genes), with an update of tissue-wide RNA-seq atlas in the Melonet-DB database. Comparison between Harukei-3 and DHL92, the first published melon genome, enabled identification of 24,758 one-to-one orthologue gene pairs, whereas others were candidates of copy number variation or presence/absence polymorphisms (PAPs). Further comparison based on 10 melon genome assemblies identified genome-wide PAPs of 415 retrotransposon Gag-like sequences. Of these, 160 showed fruit ripening-inducible expression, with 59.4% of the neighboring genes showing similar expression patterns (r > 0.8). Our results suggest that retrotransposons contributed to the modification of gene expression during diversification of melon genomes, and may affect fruit ripening-inducible gene expression. Ryoichi Yano et al. report a Nanopore-based reference genome assembly of muskmelon—a fruit known for its many cultivated varieties, including cantaloupe and honeydew—using the Japanese Harukei-3 cultivar. They identify structural genetic variation by comparing the reference to several melon genome assemblies and investigate tissue-wide gene expression patterns by RNA sequencing.

M elon (Cucumis melo L.) is one of the most economically important fruit crops in the world and is a source of vitamins, minerals, and other health-promoting substances. It is thought to have originally diversified in India and Asia and is known to exhibit very wide natural variation, especially in fruit phenotypes [1][2][3] . At least 19 horticultural subgroups and six major groups of melon have been identified. A particularly notable feature of melon is the coexistence of both climacteric (ethylene-producing and showing a burst in respiration at the onset of ripening) and non-climacteric fruit types [4][5][6] . For example, the French cultivar "Vedrantais", belonging to the subgroup var. cantalupensis, is well-known example of a climacteric melon, whereas melons of the subgroup var. inodorus (e.g., the American cultivars "Honey dew" and Spanish "Piel de Sapo") are non-climacteric. The molecular mechanism of ethylene production has been intensively studied in melon, given the importance of this hormone in regulating climacteric fruitripening traits such as shelf life, which is of considerable economic importance [7][8][9][10][11] .
The melon genome comprises 12 chromosomes, and its genome size was estimated to be~454 Mb based on the nuclear DNA content 12 . This is larger than the genomes of other cucurbit plants such as Cucumis sativus (7 chromosomes, 367 Mb) and Citrullus lanatus (11 chromosomes, 425 Mb). The first reported whole genome sequence of melon was that of the experimental line designated DHL92 13 , which was originally derived from a cross between the non-climacteric "Piel de Sapo" (subsp. melo var. inodorus) and the Korean landrace "Songwhan Charmi" (subsp. agrestis var. chinensis). A genomic DNA sequence of 417 Mb was published in the latest version of the DHL92 genome reference (CM3.6.1), of which 337 or 79. 6 Mb was the actual nucleotide sequence or ambiguous bases (e.g., NNN), respectively 14 . In addition, 29,980 protein-coding genes have been reported in the genome annotation CM4.0 14 . The DHL92 genome reference has been utilized for supporting transcriptome analyses as well as quantitative trait loci (QTL) studies of important agricultural traits, including fruit ripening, fruit morphology, and disease resistance 5,11,[15][16][17][18][19] . With decreasing costs of whole genome sequencing, several melon accessions have been sequenced and characterized using the Illumina short read next-generation sequence (NGS) platform 20,21 . However, third generation sequencing technologies (e.g., PacBio RSII/sequel and Oxford Nanopore Technology [ONT]), implemented with single molecule sequencing that can generate long reads (e.g. >10 kb) emerged as an alternative approach. Many plant genomes have been assembled and/or re-evaluated using such newer DNA sequencers [22][23][24][25][26][27][28] , including the genome of the Chinese inodorus melon cultivar Payzawat, which was sequenced using a PacBio RSII platform 29 . In addition, such approaches can provide insights into genome structural variation, as was demonstrated in Solanum lycopersicum using ONT ultra-long sequencing technology 22 .
At present, little is known about genomic structural variation, especially copy number variation (CNV) and presence/absence polymorphism (PAP) in the genomes of melon subgroups. In this current study, we assembled the whole genome sequence of the semi-climacteric Japanese cultivar "Earl's Favorite Harukei-3" (Harukei-3; var. reticulatus) by coupling ultra-long ONT sequencing (R9.4.1 + R10 flow cells) with Bionano optical mapping and Illumina mate pair sequencing. This melon shows moderate climacteric ripening behavior, although the rate of ripening is less than that of the well-studied cultivar Charentais (var. cantalupensis) 4 . ONT RNA-seq-based gene prediction coupled with other methods, such as ab initio prediction, also identified 33,829 protein-coding genes whose protein BUSCO benchmark value was 1372 (95.3%). We also expanded our tissue-wide transcriptome (RNA-seq) dataset of the Melonet-DB (https://melonet-db.dna.affrc.go.jp/ or https://gene.melonet-db.jp) by adding RNA-seq samples of ethylene-producing ripening fruit. In addition, the genomes of seven more melon accessions were sequenced and assembled by ONT at the contig level to conduct the assembly-based PAP analysis of retrotransposon Gag-like sequences. Based on a combination of comparative genomics and comprehensive transcriptome analysis, we suggest that retrotransposons played a role in the modification of gene expression as well as evolution of fruitripening-inducible gene expression during diversification of melon genomes.

Results
Genome assembly and comparative genomics of Harukei-3 melon. Melon is usually described as producing sweet fruit; however, Harukei-3 produces considerably sweeter fruit than other melon accessions if it is grown in the appropriate seasons ( Supplementary Fig. 1). Indeed, the Japanese word "Harukei" means a line suitable for growing in the spring. As a consequence of its taste and attractive appearance (Fig. 1a), Harukei-3 has been used for a long time in Japan as a standard melon to breed high-grade muskmelon. To investigate the genome structure of Harukei-3, and to obtain functional gene information for future genetic studies and breeding, we assembled its genome sequence by combining ONT ultra-long sequencing (R9.4.1 and R10), Bionano optical map, Illumina Hiseq, mate pair, and linkage map information (summarized in Supplementary Fig. 2). The ONT platform yielded more long reads than did PacBio RSII, and these ONT long reads were used to generate a contig assembly with N 50 = 8.6 Mb that was 10 times higher than that of the PacBio contig ( Fig. 1b and Supplementary Fig. 2). After scaffolding with the Bionano map and mate pair data, we obtained 80 ONT-based genomic scaffolds with an N 50 = 17.5 Mb. We also assembled the genomic scaffold based on PacBio RSII using the same procedure (N 50 = 11.4 Mb). Using PacBio-based scaffolds as a hint to modify the ONT-based scaffolds, we obtained 66 genomic scaffolds with a physical gap number of 92 and an N 50 = 18.9 Mb. Finally, 12 chromosome sequences were constructed using 28 scaffolds based on linkage map information ( Fig. 1c and Table 1; Harukei-3 ver. 1.41 pseudomolecule). The chromosomal sequence lengths without ambiguous bases (e.g., NNN) from Harukei-3 (366.7 Mb) were much longer than those from DHL92 (318.2 Mb), reflecting the lower number of physical gaps in the Harukei-3 genome sequence ( Table 1). The chloroplast genome was likely to be entirely assembled because we obtained two kinds of contigs with sequence lengths of 155 and 156 kb (Table 1). When the Harukei-3 genomic sequence was compared with the linkage map of Harukei-3 16 , we observed complete colinearity between physical positions and linkage positions (Fig. 1d), indicating that the Harukei-3 genome sequence was correctly assembled at a chromosome-scale. In contrast, when it was compared with the linkage map of other accessions or sources 13,30,31 , the physical position did not match the linkage map position of some markers ( Supplementary Fig. 3), suggesting chromosome-level structural differences between the genomes of melon accessions.
We conducted genomic alignment of Harukei-3 and DHL92 and while most of the alignment showed co-linearity, a large number of small genomic sequences were observed that might have resulted from translocation across the chromosomes between the two genomes (Fig. 1e). Additionally, the Harukei-3 assembly revealed the duplication of a large genomic block (>120 kb block repeat on chromosome 5 or 9), but this was not apparent in the DHL92 genome. We previously constructed the Harukei-3 genomic pseudomolecule based on PacBio RSII data; however, this version of the genome assembly did not resolve the large genomic block duplication ( Supplementary Fig. 4), underlining the value of the ONT ultra-long reads in resolving the genomic structure. We also conducted a genomic alignment between Harukei-3 and the recently published genome of the inodorus melon, Payzawat, which was assembled based on PacBio data 29 . This revealed a large genomic block duplication on chromosome 5 of Harukei-3 that was absent, or not assembled, in the Payzawat genome ( Supplementary Fig. 5). Other large genomic blocks were also absent in the Payzawat genome (e.g., the upper part of chromosome 8).
To compare the genomes based on gene information, we newly predicted genes in the Harukei-3 genomic sequence. We mainly used ONT RNA-seq for this purpose as ONT RNA-seq analysis can yield nearly full-length sequence mRNA molecules (summarized in Supplementary Figs. 6, 7). By combining ab initio prediction (e.g., AUGUSTUS) and short read-based prediction (e.g., Braker2) as supplementary methods, we identified 33,829 protein-coding genes in the Harukei-3 genome (33,314 nuclear genes + 515 organelle genes). Both protein BUSCO benchmark analysis 32 (ver. 3.0) and an InterProScan search 33 indicated that the Harukei-3 genome annotation represents a highly comprehensive dataset of plant protein-coding genes compared with published cucurbit genomes [26][27][28][34][35][36][37][38][39][40] (complete BUSCO = 1372 [95.3%] and InterPro ID count = 5607; Table 2 and Supplementary Data 1). In addition, we observed that our ONT-based gene prediction method gave accurate insights into exon-intron gene structure, whereas other methods (e.g., AUGUSTUS and Braker2) gave incorrect assessments of structure in some genes ( Supplementary Fig. 8). We also searched for repetitive elements in the Harukei-3 genome (e.g., simple sequence repeats and DNA/ RNA transposons). In total, 211 Mb of the Harukei-3 genomic sequence was found to correspond to repetitive elements (Supplementary Data 2).
We next analyzed the gene partner relationship between Harukei-3 and DHL92 genomes based on bidirectional BLAST-n/ p search 41 and transcript alignment analysis with the BLAST-like alignment tool 42 (Fig. 2a). Of the 33,314 nuclear genes of the Harukei-3 genome, 24,747 and 1317 genes showed one-to-one orthology or a homologous partner relationship, respectively, with DHL92 genes (Supplementary Data 3). We attached a consensus gene ID to the 24,747 orthologous genes of the Harukei-3 genome (a gene ID that starts with the "MELO3C" string and ends with the ".jh1" string) to help maintain a consistent gene nomenclature among melon genomes (Fig. 2b).
Although it is generally difficult to compare reference genomes that are generated by different sequencing technologies, 1203 genes were identified as possible candidates for CNV and PAP (Fig. 2a, Supplementary Fig. 9a  between the three genome references, the Harukei-3 reference showed slightly better alignment ratios relative to the other two genome references in the resequencing data of Harukei-3 itself, Honey dew, and Spicy ( Supplementary Fig. 10a). In contrast, the DHL92 reference showed better alignment ratios in comparison to the other two references in the case of makuwa and wild melon. Unlike the alignment ratio, the number of polymorphisms (single nucleotide polymorphisms (SNPs) and small insertions/ deletions (Indels)) that were predicted to affect protein amino acid sequence was highly variable between the Harukei-3 and DHL92 references ( Supplementary Fig. 10b). These results underlined the importance of using several genome references in the resequencing study.
Co-expression analysis of fruit-ripening-inducible genes. Since Harukei-3 fruit produce ethylene during ripening, we investigated ethylene-related gene expression in Harukei-3 and updated the tissue-wide melon RNA-seq transcriptome dataset of Yano et al. 18 by adding data derived from ethylene-emitting ripe fruit (Fig. 3a). Alignment ratios of RNA-seq reads were much higher in Harukei-3 genome reference than in DHL92 probably because the RNA-seq data were obtained from Harukei-3 itself. We identified 27,687 Harukei-3 genes with Fragments Per Kilobase of exon per Million mapped fragment (FPKM) values ≥ 0.1 ( Supplementary  Fig. 11, Supplementary Data 6 and 7). Such high-sensitive detection of gene expression level enabled high-resolution coexpression analysis, and weighted genome-wide correlation network analysis 43 (WGCNA) identified >60 co-expression clusters, including those specific to ripening fruit ( Supplementary Fig. 12). We also updated the Melonet-DB web-application tools, "Gene expression map viewer" (https://melonet-db.dna.affrc.go.jp/ap/ mvw) and "Co-expression viewer" (https://melonet-db.dna.affrc. go.jp/ap/mds) (Supplementary Figs. 13,14,and 15). In the newer dataset, up-regulation of ethylene-related genes (e.g., CmACO1, CmETR1/2, and CmNOR-NAC) in ripe fruit was observed (Fig. 3b), consistent with ethylene production by the Harukei-3 fruit. Further co-expression analysis, including not only known ethylene-related genes but also 81 NAC domain, 90 homeobox, and 42 MADS-box transcription factors, identified a coexpression cluster that was specific to fruit ripening (Fig. 4a). In the central region of this cluster, we identified an AGAMOUSlike gene, MELO3C019694.jh1, which is a homolog of Tomato AGAMOUS-LIKE 1 (TAGL1); a gene that has been shown to be involved in regulating fruit ripening 44,45 . Zhao et al. 21 also recently identified this melon gene as a candidate QTL that regulates fruit suturing. A comparison of the Harukei-3 and DHL92 genomes indicated that Harukei-3 carries a longer protein-coding transcript of MELO3C019694 relative to that of DHL92, and its expression was higher in ripe fruit than in preripe fruit ( Supplementary Fig. 16). When the genome sequences flanking MELO3C019694 were analyzed using the Harukei-3 genome reference, the upstream promoter region of MELO3C019694 was found to contain two Ty3-gypsy LTR-retroelements ( Fig. 4b; chr11 . Around this genomic region, two protein-coding sequences were identified. One of them, MELO.jh102711.1, encodes a protein sequence with retrotransposon-related protein domains such as IPR005162 (retrotransposon Gag domain), IPR013242 (retroviral aspartyl protease), and IPR000477 (reverse transcription). Together with LTR/gypsy elements, these sequences appear to function as LTR retrotransposons. According to the tissue-wide transcriptome dataset, MELO3C019694 and the neighboring retrotransposonrelated protein-coding sequences (e.g., MELO.jh102711.1) showed Table 1 Summary of assembled sequence lengths and predicted genes in the Harukei-3 genome. DNA sequence lengths of each chromosome and unanchored sequence are compared between Harukei-3 ver. 1.41 reference (left) and DHL92 CM3.6.1/4.0 reference (right). In Harukei-3, organelle genomes are assembled as one scaffold.
Harukei PAPs and the expression of retrotransposon Gag sequences.
The presence of the LTR retrotransposons in the upstream promoter sequence of MELO3C019694.jh1 (SlTAGL1 homolog) in the Harukei-3 genome prompted us to perform an enrichment analysis of CNV and PAP candidates (1,203 genes or putative protein-coding sequences) found between Harukei-3 and DHL92 genomes (Fig. 2a). For the purpose, a web-application tool, designated "GO term enrichment analysis", was developed for the Melonet-DB database (https://melonet-db.dna.affrc.go.jp/ap/got). This revealed several InterPro IDs (IPRIDs) that are usually observed in retrotransposon-related function to be enriched in the 1203 candidates (e.g., IPR005162, IPR013242, and IPR000477; Supplementary Fig. 9b). This suggested that retrotransposons have been copied or jumped across chromosomes during diversification of the melon genome. Such structural difference should be detected in the form of PAPs between genomes if genome sequences are aligned and compared at a relatively narrow range (e.g., 50-100 kb). To analyze the PAPs of retrotransposon-related sequences in such a manner, we sequenced seven more melon genomes by using the ONT R9.4.1 platform. Genome assemblies were obtained as contig datasets for Natsukei-1 (var. reticulatus), Fuyukei (var. reticulatus), and Awamidori (var. conomon) in addition to Spicy, Honey dew, Ougon-9, Awamidori, and JSS6 (wild melon). The N 50 values for the contig assemblies were more than 3.5 Mb, with a maximum value of 10.1 Mb ( Supplementary Fig. 17), indicating that they were sufficient for local genomic sequence alignment analysis. According to the Harukei-3 ver. 1.41 genome reference, there are at least 415 putative protein-coding sequences with IPR005162 (retrotransposon Gag domain). Assembly-based sequence alignment analysis between Harukei-3 and other melons indicated that there are PAPs in these retrotransposon Gag-like sequences; some were conserved between melon accessions (e.g., MELO. jh102304.1) while others were not conserved and/or present in specific melons (e.g., MELO.jh102711.1 and MELO.jh033067.1) (Fig. 5a). A hierarchical clustering analysis based on the PAP genotype dataset indicated that Natsukei-1 and Fuyukei melons are much closer to Harukei-3 compared to other accessions, which is consistent with the fact that they have the same origin (a reticulatus melon imported from the United Kingdom around 100 years ago) (Fig. 5b). By contrast, the Asian melons Awamidori, Ougon-9, and JSS6 are distant from such reticulatus melons. This result indicated that the PAP datasets obtained by assembly-based analysis reflect genomic variation between melon cultivars and accessions. Then, tissue-wide gene expression patterns were analyzed in 415 retrotransposon Gag-like sequences. Interestingly, 160 (38.6%) showed fruit-ripening-inducible expression with the highest levels in post-harvest fruit samples (Fig. 5c). One of them, MELO.jh033067.1, was present only in reticulatus melons, but it showed higher expression levels in ripening fruits than did other genes (Fig. 5a, c). As described above, the expression pattern of MELO.jh102711.1, a Gag-like sequence, is similar to that of the neighboring gene, MELO3C019694 (Fig. 4c). Thus, to investigate whether neighboring genes are co-expressed with Gaglike sequences, the correlation of gene expression was analyzed in a genome-wide manner. The result indicated that 59.4% of genes neighboring fruit-ripening-inducible Gag-like sequences also showed similar expression patterns ( Fig. 5d and Supplementary Data 8; Pearson's correlation coefficient r > 0.8). In contrast, in the case of other Gag-like sequences, the degrees of correlation between the Gag-like sequences and the neighboring genes were the same levels as the control. Heat-inducible expression of retrotransposon Gag sequences.
There is increasing evidence that the expression of plant retrotransposon-related sequences is up-regulated by abiotic and biotic stress [46][47][48] . To investigate the environmental response of Harukei-3 retrotransposon-related sequences, we performed a field transcriptome analysis of leaf samples collected weekly from plants grown in a greenhouse at the University of Tsukuba from early summer to late autumn. In Japan, midsummer is not an appropriate season for melon cultivation because the temperature inside the greenhouse sometimes exceeds 45°C. Indeed, Harukei-3 plants grown during midsummer showed severe signs of heat stress damage (Fig. 6a). In contrast, Harukei-3 plants grown during a cooler period (e.g., before midsummer or after September) had thick leaves with a dense green color, and produced sweeter melon fruit ( Fig. 6a and Supplementary Fig. 1). We generated transcriptome data corresponding to a total 75 time points and 18 independent plants (Supplementary Data 9). WGCNA clustering using this dataset indicated that the retrotransposon Gag-like sequences were co-expressed with heat shock protein genes that carry the signature of HSP20-like chaperones (Fig. 6b, c) ) were also co-expressed with these Gaglike sequences (Fig. 6c, d), indicating that melon plants grown during midsummer experienced drought stress, in addition to heat stress. Taken together, these results suggested that some retrotransposon Gag-like sequences were responsive to abiotic stress such as heat stress.

Discussion
In this study, our assembly-based genome sequence comparison clearly demonstrated that there is a substantial PAPs in retrotransposon Gag-like sequences between melon genomes. For example, one of the Gag-like sequences, MELO.jh33067.1, is present in the genome of reticulatus melons (e.g., Harukei-3, Natsukei-1, and Fuyukei) but is absent from those of some melons, such as DHL92 and Payzawat (Fig. 5a). The genetic relationship was also well reflected in the hierarchical clustering dendrogram based on the PAP genotype dataset (Fig. 5b), demonstrating that our assembly-based approach is successful to analyze genome-wide PAPs between genomes. There is growing evidence that some plant LTR retrotransposons, which are classified into the subfamilies gypsy and copia, exhibit stressinducible transcription [46][47][48] . The phytohormone ethylene has also been shown to be involved in stress-induced expression of retrotransposon in Solanum chilense 49 . If a retrotransposon is inserted into a promoter region, it may affect constituent cis-acting elements, which in turn may alter or enhance the transcriptional response of downstream genes to environmental and/or developmental signals 50 . Interestingly, most of the retrotransposon Gag-like sequences, including MELO.jh33067.1, showed the up-regulation of expression in ethylene-producing ripening fruits (Fig. 5c). Moreover, 59.4% of genes neighboring the fruit-ripening-inducible Gag-like sequences were coexpressed with the retrotransposon Gag-like sequences. Therefore, it is possible that the region of LTR retrotransposon around Gag-like sequences carry cis-acting and/or enhancer DNA elements that induce fruit-ripening-inducible gene expression. It seems also possible that ethylene production indirectly drives the transcription of Gag-like sequences. In plants, the biosynthesis of ethylene is also known to be stimulated under abiotic conditions, including heat stress 51 . Indeed, some of the Gag-like sequences were also up-regulated upon heat stress in melon (e.g., MELO. jh33067.1 in Fig. 6d). In A. thaliana, heat-inducible transcription of the retrotransposon ONSEN is regulated through small interfering RNA 52,53 , and an epigenetic regulatory mechanism may similarly be involved in the regulation of LTR retrotransposon in melon. It is possible that such mechanism is a driving force to diversify the genome not only in melon but also in other plants.
The same also may be applicable to MELO3C019694.jh1, which is a homolog of tomato AGAMOUS-like gene SlTAGL1, a regulator of fruit ripening. Recently, Zhao et al. 21 reported that MELO3C019694.jh1 gene expression levels were higher in the sutured melon Vedrantais (var. cantalupensis) than in nonsutured Piel de Sapo (var. inodorus) at 7 days after pollination 21 .
Here, we identified this gene as a possible regulator of fruit ripening as it is co-expressed with known ethylene-related fruitripening genes such as CmACO1 and CmNOR-NAC (Fig. 4a). Since Vedrantais and Piel de Sapo are well-known climacteric and non-climacteric melons, respectively, it is possible that the b a    difference in gene expression levels is associated not only with fruit suturing but also with ripening behavior (e.g., presence of ethylene production). Interestingly, the upstream promoter region of MELO3C001864.jh1 in the Harukei-3 genome also carries at least one LTR retrotransposon insertion (e.g., MELO. jh102711.1). Given that this LTR retrotransposon is absent from the genomes of some melons such as DHL92 and Payzawat (Fig. 5a), it is possible that such PAP of the retrotransposon insertion is also involved in variation of gene expression levels between melon cultivars. Additional studies, using genome editing and/or comparative transcriptomics, will clarify the detailed role of this LTR retrotransposon in fruit ripening.
In this study, we present evidence that assembly level genome comparisons can elucidate structural genomic variation, including PAPs, as well as large genomic block duplications. ONT-based assembly successfully resolved structural variation (Figs. 1e, 2a,  and 5a) and we show that ONT is also useful for genome-wide gene prediction, which is important for gene-based comparative genomics study (Fig. 2a). In particular, multiple ONT genome assemblies seem essential to analyze PAPs in a genome-wide manner (Fig. 5a). Information related to the Harukei-3 genome assembly, genome annotation, and the transcriptome dataset obtained in this study can be accessed using our updated webapplication tools in the Melonet-DB database (https://melonetdb.dna.affrc.go.jp/). Together with future updates, this database will contribute to functional genomic study of melon, especially reverse genetics study using the genome editing technique and TILLING mutant population.  JSS6 (C. agrestis), were obtained from the Genebank of the National Agriculture and Food Research Organization (NARO) in Japan. Melon plants were grown using the hydroponics method in the greenhouse of the University of Tsukuba in Japan as previously reported 18 . For genomic DNA sequencing, apexes of branched shoots were detached from plants and immediately frozen in liquid nitrogen. For tissue-wide RNA-seq study, tissues shown in Fig. 3a and Supplementary Fig. 6 were similarly obtained and frozen in liquid nitrogen. For field RNA-seq study, holepunched leaf samples were also collected in a weekly manner from July to November as shown in Fig. 6. These samples were crushed to powdery frozen samples using the Multi-beads shocker instrument (Yasui Kikai Corporation, Osaka, Japan) and stored at −80°C until use.

Methods
Genomic DNA isolation and DNA sequencing. Genomic DNA was isolated using Maxwell ® 16 Tissue DNA Purification Kit (Code No. AS1030, Promega, Wisconsin, USA). Although this kit is designed to couple with an automated DNA extraction machine, we did not use this but manually isolated genomic DNA by hand to obtain long intact DNA. Isolated genomic DNA was further subjected to the Short read Eliminator XS kit (Circulomics, Maryland, USA) to remove short DNA fragments with <5 kb. For ONT sequencing, the DNA library was prepared using the Ligation sequencing kit (Code No. SQK-LSK109, ONT, Oxford, UK) according to the manufacturer's protocol. DNA sequencing run was performed with a Nanopore Minion ® device coupled with flow cell. For the genome sequencing of Harukei-3, both R9.4.1 and R10 flow cells were used, whereas only R9.4.1 was used in the sequencing of Natsukei-1, Fuyukei, Spicy, Honey dew, Awamidori, and JSS6. ONT flow cells were repetitively used at least twice with the same DNA library. To obtain DNA sequence data, basecalling was performed using a CUDA-enabled GPU server with ONT's guppy ver. 3.3.0 software. For PacBio RSII and Illumina paired end (PE) sequencing, the outsourcing service of Macrogen Japan Co. Ltd (Kyoto, Japan). was used except the Illumina PE data of Harukei-3 that was obtained with a Hiseq-2000 sequencer at Cornel University. Illumina sequencing was performed with the 100 bp PE mode in Harukei-3 or the 150 bp PE mode in Honey dew, Spicy, Manshuu, Ougon-9, and JSS6. Illumina mate pair sequencing was performed with 5 kb and 10 kb insert libraries by using the outsourcing service of Hokkaido System Science Co. Ltd. (Sapporo, Japan).
RNA isolation and RNA-seq data acquisition. Total RNA was isolated using Maxwell ® 16 LEV Plant RNA Kit (Code No. AS1430, Promega, Wisconsin, USA) according to the manufacturer's protocol. For ONT direct RNA-seq and cDNA RNA-seq, libraries were prepared using the Direct RNA sequencing kit (Code No. SQK-RNA002, ONT, Oxford, UK) or the Direct cDNA sequencing kit (Code No. SQK-DCS108), respectively, according to the manufacturer's protocol. To obtain RNA sequence data, basecalling was performed with ONT's Guppy ver.  Fig. 3. Three different groups of queries were analyzed; all genes (control), fruit-ripening-inducible Gag-like sequences, and other Gag-like sequences (not up-regulated in ripening fruits, control). In the case of fruit-ripening-inducible Gag-like sequences, co-expression was observed for 59.4% of the neighboring genes (red asterisks).

software. For
Illumina RNA-seq with the 150 bp PE mode, the outsourcing service of Macrogen Japan Co. Ltd. was used.
Whole genome assembly. The procedures, datasets, and detailed parameters used for whole genome assembly of Harukei-3 are summarized in Supplementary Fig. 2. We used two kinds of sequence reads, ONT (R9.4.1 + R10 flow cells; 32.6 Gb) and PacBio RSII (19.5 Gb), for initial contig assembly. In the case of ONT, reads with ≥5 kb were selected and used. Contigs were separately assembled based on ONT or PacBio reads with the Canu ver. 1.8 pipeline 54 ; then errors present in the contig sequences were corrected with Pion 55 using 37 Gb of Illumina PE dataset. At this point, contig N 50 values for ONT-based and PacBio-based contig assemblies were 8.6 Mb and 0.86 Mb, respectively. Then, scaffolds were assembled based on contigs by combining methods of Bionano Irys optical map and Illumina mate pair ( Supplementary Fig. 2). For Bionano scaffolding, 86 Gb raw data were obtained using the outsourcing service of AS ONE Corp. (Osaka, Japan). They were first assembled to construct "cmap" with Irys solve ver. 3.2; then, cmap was used for both scaffolding and correction of chimeric contigs (incorrectly assembled contigs) using the same software. Scaffolding was also performed using 74 Gb of Illumina mate pair (5 kb and 10 kb insert sizes) with SSPACE ver. 3.0 56 at different "k" parameter values (from k = 160 to 80, 40, 20, 10, and 5). By using k values from larger to smaller in series, we tried to maximize the connections and minimize false positives. Chimeric scaffolds generated by the mate pair scaffolding were again corrected using Bionano cmap. At this point, scaffold N 50 values for ONT-based and PacBio-based assemblies were 17.5 Mb and 11.4 Mb, respectively. ONT-based scaffolds were further updated using PacBio-based scaffolds as a hint. In this attempt, both scaffolds were first classified into each chromosome group by using linkage map information.  Heat stress-associated co-expression cluster Fig. 6 Transcriptional responses of retrotransposon Gag-like sequences in the greenhouse from July to November. a Harukei-3 melon plants under nonstressed or heat-stressed conditions. Leaf disks were collected weekly using a hole punch. In total, 75 RNA-seq data were obtained in 18 independent Harukei-3 plants from July to November. b InterPro ID enrichment analysis in the genes that co-express with retrotransposon Gag-like sequences. Co-expression dataset is based on 75 leaf RNA-seq data obtained in the greenhouse. Heat shock protein genes are highly enriched in the genes. c Co-expression of Gag-like sequences with abiotic stress-related genes. Gag-like sequences and abiotic stress-related genes are co-expressed in one cluster (shown by red dashed circle). d Changes in gene expression in the greenhouse. Expression patterns of two retrotransposon Gag-like sequences are shown together with those of abiotic stress-related genes such as AtABF3-like and AtHSP21-like.
connect two distinct ONT-based scaffolds were identified by BLAST-n search using the following conditions: p-value ≤ 1e-150, sequence identity ≥ 99%, blast score ≥ 1000, and alignment length ≥ 5000 bp ( Supplementary Fig. 2). If both end of the PacBio-based scaffold had sequence alignments with distinct ONT-based scaffolds, it was used to connect them. Finally, the chromosome-scale pseudomolecule was constructed using 28 genomic scaffolds that were anchored and oriented by linkage map information 16,30,57 .
For the contig assembly in Natsukei-1, Fuyukei, Spicy, Honey dew, Ougon-9, Awamidori, and JSS6, ONT sequencing reads with ≥5 kb (R9.4.1 flow cell) were first subjected to the Canu ver. 1.8 pipeline. The resultant contig sequences were subjected to Racon 58 and Medaka (https://nanoporetech.github.io/medaka/) to correct erroneous bases. To further determine the candidates of chimeric contigs, ONT reads were aligned to contig sequences with minmap2 59 using the following parameter: "-a -uf -k14 -A 2 -B 4 -O 4,24 -E 2,1". Then, read depth data were obtained based on the read alignment information using the mpileup function of samtools 60 . Contigs were split at positions where the depth of ONT reads was less than four.
Gene prediction. Procedures and datasets used for genome annotation are summarized in Supplementary Figs. 6, 7. For ONT-based gene prediction, datasets of 8.2 Gb ONT direct RNA-seq and 8.8 Gb ONT cDNA RNA-seq were combined and used. They were aligned to Harukei-3 ver. 1.41 genomic sequence with Minimap2 using the following parameters "-ax splice -uf -k14." Then, transcript information with exon-intron structure was obtained by pinfish (https://github. com/nanoporetech/ont_tutorial_pinfish) using several "c" parameter values (c = 2, 3, 5, and 10; Supplementary Fig. 7). Predicted transcript sequences were obtained from the genomic sequence with gffread (https://ccb.jhu.edu/software/stringtie/ gff.shtml#gffread) based on the General Feature Format (GFF) information of Pinfish; then, they were combined with transcript sequence information of DHL92 genome annotation CM4.0. Again, Minimap2 alignment and Pinfish prediction were performed based on the combined transcript sequences to obtain the merged GFF annotation. Next, the protein-coding open reading frame (ORF) was predicted in each transcript followed by hmmsearch. The best possible ORFs were kept by selecting those with the highest sum total hmmsearch scores. If no protein domain was found in any ORF candidates, the longest ORF was kept. Transcripts were further grouped into gene units based on the position of the exon(s) and the results of the self-BLAST search (both transcript and protein sequences); transcripts were grouped if both the transcript and protein sequences had homology to each other and the positions of exon(s) on the genome were consistent between them. To further select the best-possible ORF in each gene unit, hmmsearch as well as BLAST-p search against protein sequences of 9 plant genomes were performed again (a list of the genome references used for the purpose is shown in Supplementary Fig. 7d). The ONT-based method described above predicted 31,306 protein-coding genes. A perl pipeline designated "ONT4genepredict" was developed to automatically perform the information analysis described above. It is available in Melonet-DB (https://melonet-db.dna.affrc.go.jp/ap/dnl). To evaluate the completeness of predicted genes, we used BUSCO ver. 3.0 benchmark 32 . The protein BUSCO score for the ONT-based gene dataset was 1362 (94.6%) ( Supplementary Fig. 7c). In addition to the ONT-based method, gene prediction was also performed with AUGUSTUS ver. 3.3.2 (ab initio method) 61 , Braker2 pipeline 62 , and Genome Threader 63 . For AUGUSTUS gene prediction, the parameter dataset was first trained and generated with the perl script autoAug.pl using the ONT-based annotation dataset described above. Then, genes were predicted with AUGUSTUS software using the default parameters. Braker2 was executed using the Illumina RNA-seq dataset of 45 tissue-wide samples (total 118 Gb; Fig. 3a). RNA-seq reads were first aligned to the Harukei-3 genome sequence, then the read alignment information was merged and used for Braker2 gene prediction. Transcript annotation was also obtained with StringTie 64 software based on the read alignment information. Genome Threader was executed based on the protein sequences of 10 published plant genomes (listed in Supplementary Fig. 7d). Then, EvidenceModeler (EVM, https://evidencemodeler.github.io/) was used to integrate the results of StringTie, AUGUSTUS, Braker2, and Genome Threader with the weight score setting of 10, 8, 1, and 1. EVM produced the dataset of 59,613 protein-coding genes with complete BUSCO ver. 3.0 score = 1348 (93.6%). Finally, using the EVM-based dataset as a supplementary dataset, ONT-based annotation dataset was updated to obtain 33,829 protein-coding genes (40,363 transcripts, BUSCO ver. 3.0 score = 1,372 [95.3%]). InterProScan 33 was also conducted to obtain GO and InterPro ID in each predicted protein amino acid sequence.
Comparative genomics analysis. Genomic alignment was performed with LAST 69 using the following parameters, "-e 25 -v -q 3 -j 4 -P 32 -a 1 -b 1 (for lastal)" and "-s 35 -v (for last-split)". After file format conversion via mafconvert, plot graphs comparing distinct genomes were generated with R 3.2.3 based on the result of LAST alignment. To obtain information of orthologue partners, we performed bidirectional blast searches based on both transcript and protein sequences. Because it is difficult to determine one-to-one orthologue partners when genes are duplicated in either genome, we also used the information of transcript alignment obtained with blat 42 as supporting information. When genes were paired in both bidirectional BLAST-n/p search and transcript alignment analysis, they were determined as one-to-one orthologue partners. Candidates of CNV and PAP were determined by integrating the results of bidirectional BLAST-n/p search with the information of gene position on the chromosome that could be obtained from genome annotation (GFF3 files). Information analysis described above has been automated with Perl scripts. Circos 70 was used to visualize and compare genomes.
For PAP analysis of retrotransposon Gag-like sequences, DNA alignment analysis was performed for each sequence in nine melon genomes using the Harukei-3 genome as a standard reference. First, the genomic sequence containing the Gag-like sequence and its surrounding region (approx. 50-100 kb) was obtained from the Harukei-3 genome sequence. By using this sequence as a query, a BLAST-n search was performed against the whole genome sequence data of the target melon cultivar or accession. The specific region of contig or chromosome that showed the best homology to the query sequence was identified from the BLAST search result, then the DNA sequence of this region was obtained and further used for LAST alignment. Plot graphs comparing both sequences (Harukei-3 versus the target melon) were generated with R 3.2.3 based on the result of the LAST alignment. The PAP status of the Gag-like sequence in the target melon was also calculated as a numeric value based on the alignment ratio of the corresponding genomic region. For example, an alignment ratio of 1.0 means the complete presence of the Gag-like sequence in the target melon genome, while 0 means that the sequence is absent.
InterPro ID enrichment analysis. ID enrichment analysis was performed based on Fisher's exact test using the R-exact2x2 module. To calculate q-value from p-value, Rqvalue package was used. This ID enrichment analysis is available in the "GO enrichment analysis tool" in Melonet-DB (https://melonet-db.dna.affrc.go.jp/ap/got).
Statistics and reproducibility. All statistical tests were performed using available softwares, packages, and online tools mentioned in the methods. Reproducibility can be accomplished using raw sequencing data deposited on public databases and the same command lines mentioned in the methods, where we used publicly available softwares for most of the analysis. Fisher's exact test was used for testing enriched GO or InterPro terms. Both p-value and q-value was used to indicate statistical significance. The number of RNA-seq samples used for tissuewide or leaf co-expression analyses were 45 and 75, respectively, which were determined according to the previous study 18 . Co-expression was evaluated based on the weight values calculated by R-WGCNA and pearson's correlation coefficients (n = 45 or 75).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Code availability
Web-application tools used in this study are available on Melonet-DB (https://melonetdb.dna.affrc.go.jp) or its mirror site (https://gene.melonet-db.jp/). The ONT4genepredict pipeline that can be used to predict genes based on ONT RNA-seq is also available on this web site.