ORIGINAL RESEARCH published: 31 January 2020 doi: 10.3389/fpls.2019.01815 An Improved Melon Reference Genome With Single-Molecule Sequencing Uncovers a Recent Burst of Transposable Elements With Potential Impact on Genes Raúl Castanera 1†, Valentino Ruggieri 1,2†, Marta Pujol 1,2, Jordi Garcia-Mas 1,2* and Josep M. Casacuberta 1* 1 Centre for Research in Agricultural Genomics CSIC-IRTA-UAB-UB, Campus UAB, Edifici CRAG, Barcelona, Spain, Edited by: 2 Institut de Recerca i Tecnologia Agroalimentàries (IRTA), Genomics and Biotecnology Program, Barcelona, Spain Sean Mayes, University of Nottingham, United Kingdom The published melon (Cucumis melo L.) reference genome assembly (v3.6.1) has still 41.6 Mb Reviewed by: Zhangjun Fei, (Megabases) of sequences unassigned to pseudo-chromosomes and about 57 Mb of gaps. Cornell University, Although different approaches have been undertaken to improve the melon genome assembly United States in recent years, the high percentage of repeats (~40%) and limitations due to read length have Chee Keng Teh, Sime Darby Plantation, made it difficult to resolve gaps and scaffold's misassignments to pseudomolecules, especially Malaysia in the heterochromatic regions. Taking advantage of the PacBio single- molecule real-time *Correspondence: (SMRT) sequencing technology, an improvement of the melon genome was achieved. About Jordi Garcia-Mas jordi.garcia@irta.cat 90% of the gaps were filled and the unassigned sequences were drastically reduced. A lift-over Josep M. Casacuberta of the latest annotation v4.0 allowed to re-collocate protein-coding genes belonging to the josep.casacuberta@cragenomica.es unassigned sequences to the pseudomolecules. A direct proof of the improvement reached in † These authors have contributed the new melon assembly was highlighted looking at the improved annotation of the equally to this work transposable element fraction. By screening the new assembly, we discovered many young Specialty section: (inserted less than 2Mya), polymorphic LTR-retrotransposons that were not captured in the This article was submitted to Plant previous reference genome. These elements sit mostly in the pericentromeric regions, but Breeding, a section of the journal Frontiers in Plant Science some of them are inserted in the upstream region of genes suggesting that they can have Received: 15 July 2019 regulatory potential. This improved reference genome will provide an invaluable tool for Accepted: 30 December 2019 identifying new gene or transposon variants associated with important phenotypes. Published: 31 January 2020 Keywords: long-reads, assembly, reference genome, transposable elements, melon Citation: Castanera R, Ruggieri V, Pujol M, Garcia-Mas J and Casacuberta JM (2020) An Improved Melon Reference Genome With Single-Molecule INTRODUCTION Sequencing Uncovers a Recent Burst Melon (Cucumis melo L.) is one of the most important plant crops, with a worldwide production of Transposable Elements With Potential Impact on Genes. reaching near 32 million tons in 2017 (http://www.fao.org). A high-quality reference genome Front. Plant Sci. 10:1815. assembly of melon was released in 2012 (Garcia-Mas et al., 2012). This assembly was generated doi: 10.3389/fpls.2019.01815 using 454 reads and Sanger sequencing of BAC ends and contained up to 375 Mb of sequence Frontiers in Plant Science | www.frontiersin.org 1 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome assembled into 1,594 scaffolds, with an N50 of 4.68 Mb. Since nitrogen for tissue homogenization. After isopropanol then, additional improvements have been performed. In precipitation, instead of centrifugation, the DNA was recovered particular, a high-resolution genetic map was used to anchor by fishing with a little glass hook to avoid fragmentation. We up to 98.2% of the scaffold assembly to the 2X = 24 melon added a purification step using phenol:chloroform:isoamyl chromosomes (Argyris et al., 2015), followed by an optical alcohol (25:24:1), and resuspended in Milli-Q® water. DNA mapping used to improve the orientation of the scaffolds of the integrity was evaluated by gel electrophoresis and quantified by previous assembly and accurately define the gap content Qubit 2.0. DNA was purified with AMPure® PB beads, and (Ruggieri et al., 2018). Besides the efforts done to improve the length was evaluated with the Fragment Analyzer Femto Pulse original assembly, the lastest published melon reference genome (Advanced Analytical Technologies, Inc.). DNA sequencing was (v3.6.1) (Ruggieri et al, 2018), still contains up to 19.1% of its performed using Pacific Biosciences (PacBio) RSII technology at sequence in gaps and 41.6 Mb of unassigned sequences (22,123 the Platform GENTYANE, INRA/UCA (Clermont-Ferrand, unassigned contigs out of the 42,067 contigs, grouped as Chr0). France). A total of ∼2,5 million PacBio long reads were Previous analyses on the melon genome architecture have generated, which corresponds to ∼50x coverage of the described that this species contains expanded pericentromeres estimated melon genome. arising from massive amplification of transposable elements (TEs) in the past 10 million years (Mya) (Morata et al., 2018). Genome Assembly TEs tend to accumulate in centromeric and pericentromeric The reads from the PacBio system were assembled using the regions due to the preferential insertion of some elements, hierarchical genome-assembly process 4 (HGAP4) pipeline including some retrotransposon families (Neumann et al., (Pacific Biosciences, SMRT Link Suite 6.0). The principle and 2011), but also as a consequence of the counter-selection of workflow of HGAP pipeline consists of different concatenated insertions in genic regions that are more likely to be deleterious steps, including (i) the selection of the longest reads as a seeding (Contreras et al., 2015). Due to the enriched proportion of long, sequence data set, (ii) the use of each seeding sequence as a repeated sequences such as Long Terminal Repeat (LTR) reference to recruit shorter reads and preassemble reads through retrotransposons and other TEs, plant centromeres, and a consensus procedure, (iii) the assembly of the preassembled pericentromeres are difficult to assemble and often contain reads, (iv) the refinement/polishing using all initial read data to multiple gaps. Such difficulty arises from the limitation of the generate the final consensus (Chin et al., 2013). In the pre- short-read sequencing to distinguish between near-identical assembly step, the raw reads were filtered using default settings repeated regions. However, this limitation will also make it with read quality (rq) of ≥ 0.65. Then, the assembly step was difficult to correctly assemble TEs sitting in gene-rich regions. performed using FALCON in the HGAP4 tool with seed As a consequence, short-read-based assemblies may contain an coverage set to 50, “aggressive” option turned on, and underestimated transposon content, with elements missing in minimum accuracy set to 65. The ARROW algorithm was used the pericentromeric regions, but also in the proximity of genes, to polish the genome assembly with default parameters. and potentially impacting on gene regulation or coding capacity. Third generation sequencing offers a great opportunity to Reference-Guided Contig Ordering, improve short-read-based assemblies such as the melon Orientation, and Genome Quality reference genome due to the longer read length, the low Assessment systematic bias, high consensus read accuracy, and improved The contigs produced by the assembly were ordered and oriented assembling algorithms. In the recent years, many studies have based on the latest melon assembly (v3.6.1) with the RaGOO taken advantage of these technologies for improving draft tool, which uses a reference-guided process (Alonge et al., 2019). genomes or generating chromosome-level assemblies (Jiao In order to improve the mappability of PacBio contigs, a et al., 2017b; Zhang et al., 2019). We describe here a new polishing step was previously performed on the v3.6.1 reference assembly for the cultivated melon (v4.0). This new assembly using raw PacBio reads. With this aim, the ARROW version benefited from ~50-fold PacBio long reads coupled with pipeline in the SMRT Link suite (resequencing pipeline) was a 20-fold Illumina short-reads data, which allowed to improve used with default parameters, just superimposing the minimum the characterization and accuracy of several regions of the number of reads to call variant ≥ 15. genome, particularly repetitive regions and centromeric areas. RaGOO is an open-source tool, implemented as a python command-line utility, which internally invokes Minimap2 (Li, 2018). Default parameters were used with k-mer size and METHODS window size both set to 19 bp. Any alignment shorter than 1 kbp in length was removed. As reported by the author, to cluster DNA Extraction and Sequencing contigs, the tool assigns each contig to the reference chromosome Genomic DNA was extracted from the double-haploid line which it covers the most. Subsequently, for each pseudomolecule DHL92, the same line sequenced to obtain the previous group, the contigs in that group are ordered and oriented relative version of the melon genome, v3.6.1 (Ruggieri et al., 2018), as to each other by examining the longest (primary) alignment. described by (Doyle, 1991) with minor modifications. Three Ordering is then achieved by sorting these primary alignments. grams of young leaves were harvested and frozen in liquid To produce pseudomolecules, the contigs are concatenated, with Frontiers in Plant Science | www.frontiersin.org 2 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome padding of 1,000 “N” characters placed between contigs. Finally, using the Arabidopsis mutation rate of 7x10-9 nucleotides per site the new consensus sequences were polished with 20× Illumina per year (Ossowski et al., 2010). paired-end reads (2 × 150 bp). Reads were aligned to the assembly using BWA-MEM (Li and Durbin, 2009). Sequence Identification of Polymorphic LTR- error correction was performed with the Pilon pipeline (Walker Retrotransposons et al., 2014). The completeness of the final assembly was Resequencing short-read data from six melon varieties (CV, evaluated with BUSCO (version 3) (Simão et al., 2015) using IRAK, PI 161375, Trigonous, Calcuta, and Vedrantais) were the conserved plant genes as database (Eudicotyledons odb10*). previously available (Garcia-Mas et al., 2012; Sanseverino et al., Comparative analysis and synteny between v3.6.1 and v4.0 2015). Trimming and adapter removal was performed with assemblies were performed using MAUVE (Darling et al., AdapterRemoval (Lindgreen, 2012). Clean reads were mapped 2004) and SyMAP v4.2 (Soderlund et al., 2011). to the v4.0 assembly using BWA-MEM (Li and Durbin, 2009). PINDEL (Ye et al., 2009) was run on mapping files to identify Genome Annotation deletions in re-sequenced varieties, using the following The genome annotation was performed by transferring through parameters: Minimun mapping quality = 35, minimum a liftover process the latest published gene models (Ruggieri number of supporting reads for calling a deletion = 5. A et al., 2018) to the new PacBio-based genome assembly using polymorphic insertion was scored when the deletion and Maker v2 program (Campbell et al., 2014). The parameters used reference element displayed a reciprocal intersect of 90% of in the configuration file were the following: est_forward = 1, the length. est2genome = 1, split_hit = 20000, min_intron = 20, single_exon = 1, single_length = 149, correct_est_fusion = 1. In case of a gene mapping on different positions of the genome, only the match with the highest Maker score was retained. Gene RESULTS ontology (GO) enrichment analysis was carried out using Genome Assembly Workflow GOATOOLS (Klopfenstein et al., 2018) The approach followed relies on a combination of different pipelines and resources as highlighted in Supplementary Annotation of Transposable Elements Figure S1. The first step took as input four PacBio runs, which Transposable elements were detected in the new genome yielded about 21 Gigabases of sequence (corresponding assembly using TEdenovo pipeline from the REPET package approximately to a 50x melon genome coverage) with an (Flutre et al., 2011), excluding the structural search. Consensus average read length of 8 kbp and an N50 of about 15 kbp sequences representing each TE family were classified into TE (Supplementary Table S1). At the pre-assembly stage, 1,499,406 orders using PASTEC (Hoede et al., 2014) and annotation of TE seed reads were selected with an average length of 12 Kbp. The copies was carried out by TEannot using two iterations. After the seed reads were used to produce about 1,469,624 pre-assembled first TEannot run, only consensus sequences that had a full- reads (Supplementary Table S1). The final consensus assembly length match in the genome were retained. A second iteration of yielded 1,178 contigs with a N50 of 714 kbp and a total genome TEannot using these consensuses was used to obtain the final size of 357.64 Mbp. The final mean coverage obtained and the annotation. We used blastx (Repbase peptide database (Bao et al., realigned subread concordance are illustrated in Supplementary 2015), cut off e-value = e-5) to identify TIR-TE copies that Figure S2. retained coding potential. Pseudomolecule Construction Specific Annotation of LTR- Following a reference-guided process, the contigs produced by Retrotransposons the assembly were ordered and oriented based on the latest LTR-retrotransposon candidates were detected by a structural melon assembly (v3.6.1) (Supplementary Figure S1). In order to approach using LTRharvest (Ellinghaus et al., 2008). Every improve the mappability of the PacBio contigs on the genome, a element was translated to the six possible frames and scanned polishing/correction step using the complete set of PacBio reads for LTR-retrotransposon-specific domains using hmmscan was undertaken. A total of 648,906 variants (271,290 deletions, (Eddy, 2011). Elements without coding potential were filtered 293,163 insertions and 84,453 substitutions) on the published out, and the remaining elements were classified into Copia and v3.6.1 genome were corrected, leading to 1% improvement in Gypsy superfamilies based on the order of the internal coding mapping of the PacBio contigs. During pseudomolecule domains, as defined by (Xiong and Eickbush, 1990). Elements construction, we could assign 96% of the contigs to the 12 lacking one or more domains were tagged as “unclassified”. melon chromosomes, leaving only 44 unassigned short contigs (average length of 7.2 kbp). Insertion Age of LTR-Retrotransposons The LTR regions of every coding element were extracted and Further Polishing of Pseudomolecules aligned with MUSCLE (Edgar, 2004). Kimura 2P distance of The last step of the workflow was aimed to correct/polish the every aligned LTR pair was calculated and used to estimate PacBio assembly using 20-fold Illumina reads from a previous insertion age as previously reported (SanMiguel et al., 1998), study (Sanseverino et al., 2015). A total of 169,279 variants Frontiers in Plant Science | www.frontiersin.org 3 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome FIGURE 1 | Comparison of the chromosomes length (ungapped) between the v3.6.1 and the v4.0 genome assemblies. (117,320 insertions, 29,353 deletions, 22,606 substitutions), assembly has captured most of the gene information. To representing less than 0.1% of the total genome size in length, maintain gene models and names, the current annotation was were detected and included in the final genome v4.0. After error transferred to the v4.0 PacBio assembly through a liftover correction using short sequence reads, the total size of melon process. We successfully moved 28,299 out of 29,980 gene pseudomolecules is 358 Mb. The new reference assembly models to the new genome assembly. The 5% (1,618) of contains 1,169 artificial gaps (strings of 1,000 Ns) and has a transcripts that did not pass the MAKER's thresholds to define much higher contiguity than the previously published short-read a proper gene model mainly consist of proteins with unknown genome assembly DHL92 v3.6.1 (contig N50 improved from function (651), transposons (79), and girdin-like proteins (39). 26.1 kb to 714 kb; contig number improved from 42,067 to In terms of distribution, about 22% of them (374) were from the 1,178). Figure 1 shows an improvement of the v4.0 genome unassembled contigs in Chr0. This failure could be due to the fact assembly in terms of length increase of each chromosome and that part of these genes, especially those with unknown reduction of unassigned contigs in Chr0 reduction when functions, may represent false or partial gene models in the compared with v3.6.1 assembly. These results highlight the previous genome annotation. The re-arrangements of some increase of the pseudomolecules sizes of about 40 Mbp (about contigs in the new assembly may also be in part responsible of 20 Mbp already present in Chr0 and 20 Mbp of completely new this discrepancy. A complete list of these genes is provided in sequence) in this new assembly, which corresponds Supplementary Table S3. approximately to 11% of the total genome length. Supplementary Table S2 reports the anchoring of 21,283 Assembly v4.0 Captures a Larger Fraction unassigned contigs of the v3.6.1 Chr0 (96.2%) on the new of Repetitive Elements PacBio melon assembly. A synteny analysis of v3.6.1 and v4.0 We used the TEdenovo pipeline from the REPET package to assemblies showed a high degree of correspondence across all identify TE sequences from the melon v4.0 assembly and to build chromosomes, with short re-oriented or reordered blocks on all TE consensuses. Subsequently, two iterations of TEannot were chromosomes except in Chr3, Chr8, and Chr11 (Supplementary run to annotate TE sequences. Transposons cover 45.2% of the Figures S3 and S4). In terms of block relocation among new genome assembly (excluding unclassified sequences), in chromosomes noteworthy changes were detected between comparison to 35.7% found in the v3.6.1 (Morata et al., 2018). Chr02, Chr11, and Chr12 of v4.0 assembly and Chr05, Chr06, Similarly to what was found in the v3.6.1 genome assembly, LTR- Chr08, and Chr10 of the v3.6.1 assembly, respectively retrotransposons represented the largest fraction of TEs in the (Supplementary Figure S3). v4.0 assembly, followed by Terminal Inverted Repeats (TIRs) In order to assess the level of completeness of the new containing TEs (TIR-TEs) (Table 1). The number of LTR- assembly with respect to the gene content, we performed a retrotransposons is higher for v4.0 but the genome fraction BUSCO analysis. We obtained 94.8% of complete and 1.7% of that LTR-retrotransposons account for in the two assemblies is fragmented BUSCOs at the genome level and a 91.1% of similar. On the contrary, we observed a drastic increase in the complete and 2.4% fragmented BUSCOs at the gene model amount of annotated TIR-TEs and the genome fraction they level. The observed values are comparable to the ones reported account for in v4.0 as compared with v3.6.1 (14.97% and 7.11% in the v3.6.1 genome assembly, suggesting that the previous respectively). We tested for TIR-TE copies that retained coding Frontiers in Plant Science | www.frontiersin.org 4 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome TABLE 1 | Comparison of the transposable element (TE) annotation based on supports this result, as the biggest differences between the two the v3.6.1 and v4.0 assemblies. annotations are found for sizes between 100–500 bp, which are TE order Acronym v3.6.1 v4.0 compatible with the length of MITEs and partial TE copies. These differences can be attributed in part to different annotation Copies Genome Copies Genome thresholds. Indeed, the peak found at 100 bp in v4.0 absent in fraction (%) fraction (%) v3.6.1 reflects a difference in the annotation approach (minimum LTR RLX 74,161 23.44 136,761 23.81 annotation size = 200bp in v3.6.1). However, v4.0 contains more LINE RIX 11,913 2.64 15,067 1.80 TIR-TEs elements of all sizes, and in particular of elements with a SINE RSX 391 0.04 746 0.04 size shorter than 1,000 nt that probably represent truncated TIR- DIRS RYX 4,212 1.65 17,890 4.39 TIR DTX 21,383 7.11 92,819 14.97 TEs and MITEs, which could be the result of a more complete Helitron DHX 1,699 0.3 5,637 0.45 representation of repetitive sequences in the assembly. Others 912 0.48 823 0.07 TOTAL 35.66 45.53 The TE orders are referred to as follows: LTR-retrotransposons (LTR), Long Interspersed Assembly v4.0 Contains Many Young LTR- nuclear elements (LINE), short interspersed nuclear elements (SINE), DIRS retrotransposons (DIRS), TIR-TEs (TIR), and Helitrons. Retrotransposons Missing In v3.6.1 Among the different classes of TEs that populate plant genomes, young LTR-retrotransposons are the most difficult potential and found that in both cases the vast majority of the to assemble due to their length and high similarity between annotated elements were non-coding (89.1% in v3.6.1 and 96.1% copies. LTR-retrotransposons are frequently abundant and in v4.0). The v4.0 assembly has 851 more coding TIR-TEs as show a high level of polymorphism in varieties and compared with v3.6.1. However, the main difference between the individuals, which make them important targets of study. In TIR-TE fraction of both assemblies is explained by the order to annotate these elements in the v4.0 assembly and differential amount of non-coding elements. The size compare the LTR-retrotransposon content with that of v3.6.1, distribution of TIR-TEs (Supplementary Figure S5) also we used a structural and homology-based approach to identify and to date LTR-retrotransposon insertions. Using this approach, we annotated 1,320 full-length elements more in TABLE 2 | Annotation of full-length LTR-retrotransposons. Number of full-length v4.0 than in v3.6.1, which represents an increase of 40% retrotransposon copies belonging to Gypsy, Copia, and unclassified superfamilies (Table 2). An important fraction of the new LTR- in the published v3.6.1 and the v4.0 genome assemblies. retrotransposons belongs to the Gypsy superfamily, but the Superfamily v3.6.1 v4.0 v4.0 assembly also contains more Copia LTR-retrotransposons than the v3.6.1. We dated the insertion of all full-length LTR- Gypsy 815 1,526 Copia 1,067 1,427 retrotransposons, and the results showed that the vast majority Unclassified 1,358 1,607 of newly assembled elements in v4.0 are very young, with TOTAL 3,240 4,560 estimated insertion times from 0 to 2 Mya (Figure 2). The FIGURE 2 | Distribution of insertion age of Gypsy, Copia, and unclassified LTR-retrotransposons annotated on the genome assemblies v3.6.1 and v4.0. Frontiers in Plant Science | www.frontiersin.org 5 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome distribution of v4.0-specific LTR-retrotransposons followed the therefore may be affecting the promoters of such genes. In general distribution profile of all the annotated TEs, with an addition, almost 60% of these elements (69) were predicted to accumulation along the centromeres and pericentromeres in all be polymorphic in the six varieties analyzed (Supplementary the chromosomes (Figure 3), their abundance decreasing in Table S5). An example is the polymorphic Gypsy LTR- gene-rich regions. retrotransposon inserted into the promoter of an AGAMOUS MADS box transcription factor (MELO3C000260, Figure 4). Young LTR-Retrotransposons Are Highly This element is young (0.2 Mya), and predicted to be absent in 4 Polymorphic and Have a Potential Impact out of the 6 re-sequenced varieties. This Gypsy element could not on Genes be properly assembled in v3.6.1, which shows several gaps in the In order to determine to what extent the new assembled elements corresponding region upstream the AGAMOUS gene. In missing in v3.6.1 had a potential impact on genes, we analyzed in addition, in the v3.6.1 assembly the gene was located in the detail all the young, full-length LTR-retrotransposons (0-2 Mya). artificial Chr0, which contained the unassembled contigs, V3.6.1 contains 443 of these elements, whereas v4.0 contains up whereas in v4.0 assembly we could place it in Chr11 at the to 1,523. Using resequencing data from six varieties and the new position 23,043,868 bp-23,044,427 bp. A manual inspection of v4.0 assembled genome as a reference, we were able to determine this region allowed the correction of the AGAMOUS MADS box the level of polymorphism of these elements. More than half transcription factor gene by combining both MELO3C0002360 (777) of these young elements were predicted to be absent in at and MELO3C019694 (Supplementary Figure S6). Other least one of the six varieties. The newly discovered young LRT- examples of genes carrying a newly assembled LTR- retrotransposons in the genome assembly v4.0 were further retrotransposon absent in v3.6.1 are a TMV resistance protein studied for their potential impact on genes. We found that 116 N-like (MELO3C021852.2) and a UV radiation resistance- out of the 1,523 were located in the close upstream regions of associated protein (MELO3C020442), among others (Figure 4, annotated genes (< 1,000 bp, Supplementary Table S4), and Supplementary Table S4). A gene ontology enrichment analysis FIGURE 3 | Distribution of transposable elements (TEs) and genes across v4.0 pseudomolecules. In green, density of REPET features per window (6,000 windows in total). In red, density of full-length Long Terminal Repeat (LTR)-retrotransposons annotated in the v4.0 that were absent in the v3.6.1 assembly. In orange, density of polymorphic LTR-retrotransposons with insertion time below 2 Mya. In purple, gene density. Frontiers in Plant Science | www.frontiersin.org 6 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome FIGURE 4 | Example of LTR-retrotransposon insertion (red boxes) in the proximal upstream region of genes (blue boxes) annotated in the V4.0 assembly and that corresponded to a gap in the v3.6.1 assembly. found no enriched terms within the functional annotation of can complicate the interpretation of GWAS. All these limitations the 116 genes. of incomplete assemblies for genotype to phenotype studies have been previously highlighted (International Wheat Genome Sequencing Consortium (IWGSC) et al., 2018; Benevenuto et al., 2019). Here, we combined PacBio (50x) with Illumina DISCUSSION (20x) reads to improve the genome assembly v3.6.1 of the melon An Improved Melon Reference Genome reference genome DHL92. The use of 2nd generation Illumina sequencing technology to correct PacBio long reads is reported to Assembly Produced Using Long-Read be an efficient and cost effective way to improve a genome Sequencing assembly (Mahmoud et al., 2019). This integrated workflow A high-quality and accurate reference genome represents a produced an improvement of the genome assembly both in relevant resource for basic and applied research including terms of new sequence gained (20 Mbp) and inclusion of functional genetics, comparative genomics, and population previous unassigned contigs (20 Mb). In addition, short blocks genetics (Kingan et al., 2019). Indeed, many reference genomes were reoriented or reordered within and across chromosomes. for crop plants have been generated over the past decade, even The structure of the genome is therefore improved in assembly though most of them are often fragmented and missing complex version v4.0 presented here. The assessment of genome repeat regions (Jiao et al., 2017b). Melon is a widely cultivated completeness and sequence accuracy of the v4.0 assembly was crop in the world, and its reference genome was first published in performed using a set of “Eudicotyledon” conserved genes. This 2012 (Garcia-Mas et al., 2012). This reference genome has been analysis indicated that the main assembly improvements improved over the time (Argyris et al., 2015; Ruggieri et al., occurred in non-genic regions, in line with what has already 2018), and consists of 42,067 small contigs, assembled in 13 been reported for other genomes (Jiao et al., 2017a). The number scaffolds. Some of the contigs were still arbitrarily ordered and of gaps was reduced from 44,650 in the genome version v3.6.1 to oriented, which complicated the analysis of some individual loci. 1,169 in v4.0 (Supplementary Figure S7). The remaining gaps In addition, the last published version of the assembly, version probably result from the presence of very complex regions in the v3.6.1, also contains a high number of short gaps, frequently genome that will need further efforts to be solved. found in intergenic regions and often close to genes. These drawbacks are a limitation for genotype to phenotype analyses, as gaps may contain sequence variability that cannot be used for V4.0 Assembly Uncovers a Burst of Young GWAS or fine-mapping studies, and can also contain candidate LTR-Retrotransposons genes that cannot be associated to the trait. On the other hand, The new v4.0 assembly of the melon genome contains up to 10% finding significant SNPs scattered across unassigned scaffolds more TE content than v3.6.1, a difference that can be explained Frontiers in Plant Science | www.frontiersin.org 7 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome mainly by a better capture of coding and non-coding TIR-TEs, as retrotransposons in gene-rich regions. The detection of recent well as a better identification of young LTR-retrotransposons. LTR-retrotransposon insertions at close proximity of genes (< 1 Non-coding TIR-TEs such as MITEs have been described to be kb away) indicated potential to alter or regulate gene expression. involved in gene regulation through the amplification of In this study, 60% of these LTR-retrotransposons were first transcription factor binding sites (TFBS) (Hénaff et al., 2014; found to be polymorphic in the six melon varieties, providing Morata et al., 2018). Thus, our new dataset represents a probable association with phenotypic variation in melon species. significant improvement that can be used to assess the Further studies should be addressed to demonstrate this functional impact of these elements with a much better hypothesis. One of these young Gypsy LTR-transposons is precision. Besides the importance of TIR-TEs, LTR- inserted into the promoter of the AGAMOUS MADS box retrotransposons are the most interesting TEs due to their high transcription factor MELO3C019694, which was miss- abundance and their potential impact on genes. In this v4.0 annotated in assembly v3.6.1 (Supplementary Figure S6). assembly, the LTR-retrotransposon content is similar to that of Recently, MELO3C019694 has been suggested as the candidate v3.6.1 in percentage of genome fraction. However, we found a gene for the presence of sutures trait after performing GWAS large difference in the content of full-length and young elements. and bi-parental mapping experiments (Zhao et al., 2019), and the Full LTR-retrotransposons are difficult to annotate with orthologous SHP1 and SHP2 in Arabidopsis regulates pod approaches that use genome self-comparison followed by dehiscence in this plant (Liljegren et al., 2000). The insertion RepeatMasker annotation (i.e., as REPET does). This approach of a Gypsy element in the promoter of MELO3C019694 will have can be effectively used to identify truncated and degenerated to be tested in a wide collection of non- and sutured accessions. copies, but often leads to the fragmentation of long intact elements. To overcome this problem, a structural detection (LTRharvest) followed by a homology-based approach was used to identify full-length elements with coding potential in CONCLUSION both v3.6.1 and v4.0 assemblies. Using the same annotation We present here a new assembly of the melon genome, based on pipeline, we were able to identify up to 40% more full-length a combination of PacBio and Illumina sequencing, with an LTR-retrotransposons in the new v4.0 assembly, an important improved sequence content and continuity with respect to the fraction of which are located in centromeric and pericentromeric previous published assembly version. The v4.0 genome assembly regions. It is well known that Gypsy elements tend to integrate in enables identification of important recent LTR-retrotransposon such regions, which are highly repetitive and difficult to insertions at genes and their polymorphism among melon assemble. In this sense, this result evidences that v4.0 assembly varieties. These insertions may affect the coding capacity or the captures a much larger fraction of the pericentromeres than expression of melon genes and may be linked to phenotypic v3.6.1 due to the improved assembly of LTR-retrotransposons, variability in agronomic traits, as for example, the presence of especially the younger ones. Our results evidence that a recent sutures in the fruit. (less than 2 Mya) burst of LTR-retrotransposons occurred in the melon genome, which was overlooked in previous analysis due to the incompleteness of the reference assembly. Based on our results, we warn that comparisons of LTR-retrotransposon DATA AVAILABILITY STATEMENT content and distribution between genome assemblies of very The raw sequencing data and the assembly are available at the different quality could be strongly biased and should be European Nucleotide Archive, ENA PRJEB34181. The assembly carefully discussed. and the gene and TE annotations are available at the Melonomics database (www.melonomics.net). Impact of Young LTR-Retrotransposons on Genes The genome of melon has been described to have recently expanded pericentromeric regions resulting from a massive TE AUTHOR CONTRIBUTIONS amplification (Garcia-Mas et al., 2012). The number of young LTR-retrotransposons found in these regions supports the JG-M and JC conceived the project. MP obtained the DNA. VR hypothesis that the pericentromeric expansion of melon and RC obtained and analyzed the data. VR, RC, JG-M, and JC occurred after the split with cucumber, which was dated about drafted the manuscript. All authors revised and approved 10 Mya (Sebastian et al., 2010). Here, we have annotated a higher the manuscript. amount of young LTR-retrotransposons (< 2 Mya) located in pericentromeric regions providing an additional support to the hypothesis of the expansion of these regions through the FUNDING accumulation of LTR-retrotransposon insertions. In addition to the important number of previously unassembled LTR- This work was supported by the Spanish Ministry of Economy retrotransposons sitting in the pericentromeric regions, the and Competitiveness grant AGL2015-64625-C2-1-R to JG-M v4.0 assembly also contains a high number of new LTR- and AGL2016-78992-R to JC, as well as Severo Ochoa Frontiers in Plant Science | www.frontiersin.org 8 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome Programme for Centres of Excellence in R&D 2016-2010 (SEV- SUPPLEMENTARY MATERIAL 2015-0533) and the CERCA Programme/Generalitat de Catalunya to both groups. RC was recipient of a Juan de la The Supplementary Material for this article can be found online Cierva Postdoctoral fellowship from the Spanish Ministerio de at: https://www.frontiersin.org/articles/10.3389/fpls.2019.01815/ Economia y Competitividad. full#supplementary-material REFERENCES assemblies of three plant species using optical mapping and chromosome conformation capture data. Genome Res. 27, 778–786. doi: 10.1101/ Alonge, M., Soyk, S., Ramakrishnan, S., Wang, X., Goodwin, S., Sedlazeck, F. J., gr.213652.116 et al. (2019). Fast and accurate reference-guided scaffolding of draft genomes. Jiao, Y., Peluso, P., Shi, J., Liang, T., Stitzer, M. C., Wang, B., et al. (2017b). BioRxiv Genome Biol. 20, 224. doi: 10.1186/s13059-019-1829-6 Improved maize reference genome with single-molecule technologies. Nature Argyris, J. M., Ruiz-Herrera, A., Madriz-Masis, P., Sanseverino, W., Morata, J., 546, 524–527. doi: 10.1038/nature22971 Pujol, M., et al. (2015). Use of targeted SNP selection for an improved Kingan, S. B., Heaton, H., Cudini, J., Lambert, C. C., Baybayan, P., Galvin, B. D., anchoring of the melon (Cucumis melo L.) scaffold genome assembly. BMC et al. (2019). A high-quality de novo genome assembly from a single Genomics 16, 4. doi: 10.1186/s12864-014-1196-3 mosquito using PacBio sequencing. Genes (Basel) 10 (1), 62. doi: 10.3390/ Bao, W., Kojima, K. K., and Kohany, O. (2015). Repbase Update, a database of genes10010062 repetitive elements in eukaryotic genomes. Mob. DNA 6, 11. doi: 10.1186/ Klopfenstein, D. V., Zhang, L., Pedersen, B. S., Ramírez, F., Warwick Vesztrocy, A., s13100-015-0041-9 Naldi, A., et al. (2018). GOATOOLS: A Python library for Gene Ontology Benevenuto, J., Ferrão, L. F. V., Amadeu, R. R., and Munoz, P. (2019). How can a analyses. Sci. Rep. 8, 10872. doi: 10.1038/s41598-018-28948-z high-quality genome assembly help plant breeders? Gigascience 8 (6), giz068. Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with doi: 10.1093/gigascience/giz068 Burrows-Wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/ Campbell, M. S., Holt, C., Moore, B., and Yandell, M. (2014). Genome annotation bioinformatics/btp324 and curation using MAKER and MAKER-P. Curr. Protoc. Bioinf. 48, 4.11.1– Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences. 4.1139. doi: 10.1002/0471250953.bi0411s48 Bioinformatics 34, 3094–3100. doi: 10.1093/bioinformatics/bty191 Chin, C.-S., Alexander, D. H., Marks, P., Klammer, A. A., Drake, J., Heiner, C., Liljegren, S. J., Ditta, G. S., Eshed, Y., Savidge, B., Bowman, J. L., and Yanofsky, M. et al. (2013). Nonhybrid, finished microbial genome assemblies from long-read F. (2000). SHATTERPROOF MADS-box genes control seed dispersal in SMRT sequencing data. Nat. Methods 10, 563–569. doi: 10.1038/nmeth.2474 Arabidopsis. Nature 404, 766–770. doi: 10.1038/35008089 Contreras, B., Vives, C., Castells, R., and Casacuberta, J. M. (2015). “The impact of Lindgreen, S. (2012). AdapterRemoval: easy cleaning of next-generation transposable elements in the evolution of plant genomes: from selfish elements sequencing reads. BMC Res. Notes 5, 337. doi: 10.1186/1756-0500-5-337 to key players,” in Evolutionary Biology: Biodiversification from Genotype to Mahmoud, M., Zywicki, M., Twardowski, T., and Karlowski, W. M. (2019). Phenotype. Ed. P. Pontarotti (Cham: Springer International Publishing), 93– Efficiency of PacBio long read correction by 2nd generation Illumina 105. doi: 10.1007/978-3-319-19932-0_6 sequencing. Genomics 111, 43–49. doi: 10.1016/j.ygeno.2017.12.011 Darling, A. C. E., Mau, B., Blattner, F. R., and Perna, N. T. (2004). Mauve: multiple Morata, J., Tormo, M., Alexiou, K. G., Vives, C., Ramos-Onsins, S. E., Garcia-Mas, alignment of conserved genomic sequence with rearrangements. Genome Res. J., et al. (2018). The evolutionary consequences of transposon-related 14, 1394–1403. doi: 10.1101/gr.2289704 pericentromer expansion in melon. Genome Biol. Evol. 10, 1584–1595. doi: Doyle, J. (1991). “DNA protocols for plants,” in Molecular techniques in taxonomy. 10.1093/gbe/evy115 Eds. G. M. Hewitt, A. W. B. Johnston and J. P. W. Young (Berlin, Heidelberg: Neumann, P., Navrátilová, A., Koblížková, A., Kejnovský, E., Hřibová, E., Hobza, Springer Berlin Heidelberg), 283–293. doi: 10.1007/978-3-642-83962-7_18 R., et al. (2011). Plant centromeric retrotransposons: a structural and Eddy, S. R. (2011). Accelerated profile HMM searches. PloS Comput. Biol. 7, cytogenetic perspective. Mob. DNA 2, 4. doi: 10.1186/1759-8753-2-4 e1002195. doi: 10.1371/journal.pcbi.1002195 Ossowski, S., Schneeberger, K., Lucas-Lledó, J. I., Warthmann, N., Clark, R. M., Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy Shaw, R. G., et al. (2010). The rate and molecular spectrum of spontaneous and high throughput. Nucleic Acids Res. 32, 1792–1797. doi: 10.1093/nar/ mutations in Arabidopsis thaliana. Science 327, 92–94. doi: 10.1126/ gkh340 science.1180677 Ellinghaus, D., Kurtz, S., and Willhoeft, U. (2008). LTRharvest, an efficient and Ruggieri, V., Alexiou, K. G., Morata, J., Argyris, J., Pujol, M., Yano, R., et al. (2018). flexible software for de novo detection of LTR retrotransposons. BMC Bioinf. 9, An improved assembly and annotation of the melon (Cucumis melo L.) 18. doi: 10.1186/1471-2105-9-18 reference genome. Sci. Rep. 8, 8088. doi: 10.1038/s41598-018-26416-2 Flutre, T., Duprat, E., Feuillet, C., and Quesneville, H. (2011). Considering SanMiguel, P., Gaut, B. S., Tikhonov, A., Nakajima, Y., and Bennetzen, J. L. (1998). transposable element diversification in de novo annotation approaches. PloS The paleontology of intergene retrotransposons of maize. Nat. Genet. 20, 43– One 6, e16526. doi: 10.1371/journal.pone.0016526 45. doi: 10.1038/1695 Garcia-Mas, J., Benjak, A., Sanseverino, W., Bourgeois, M., Mir, G., González, V. Sanseverino, W., Hénaff, E., Vives, C., Pinosio, S., Burgos-Paz, W., Morgante, M., M., et al. (2012). The genome of melon (Cucumis melo L.). Proc. Natl. Acad. et al. (2015). Transposon insertions, structural variations, and snps contribute Sci. U.S.A. 109, 11872–11877. doi: 10.1073/pnas.1205415109 to the evolution of the melon genome. Mol. Biol. Evol. 32, 2760–2774. doi: Hénaff, E., Vives, C., Desvoyes, B., Chaurasia, A., Payet, J., Gutierrez, C., et al. 10.1093/molbev/msv152 (2014). Extensive amplification of the E2F transcription factor binding sites by Sebastian, P., Schaefer, H., Telford, I. R. H., and Renner, S. S. (2010). Cucumber transposons during evolution of Brassica species. Plant J. 77, 852–862. doi: (Cucumis sativus) and melon (C. melo) have numerous wild relatives in Asia 10.1111/tpj.12434 and Australia, and the sister species of melon is from Australia. Proc. Natl. Hoede, C., Arnoux, S., Moisset, M., Chaumier, T., Inizan, O., Jamilloux, V., et al. Acad. Sci. U.S.A. 107, 14269–14273. doi: 10.1073/pnas.1005338107 (2014). PASTEC: an automatic transposable element classification tool. PloS Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V., and Zdobnov, E. One 9, e91929. doi: 10.1371/journal.pone.0091929 M. (2015). BUSCO: assessing genome assembly and annotation completeness International Wheat Genome Sequencing Consortium (IWGSC) and IWGSC with single-copy orthologs. Bioinformatics 31, 3210–3212. doi: 10.1093/ RefSeq principal investigators:, Appels, R., Eversole, K., Feuillet, C., Keller, B., bioinformatics/btv351 et al. (2018). Shifting the limits in wheat research and breeding using a fully Soderlund, C., Bomhoff, M., and Nelson, W. M. (2011). SyMAP v3.4: a turnkey annotated reference genome. Science 361, eaar7191. doi: 10.1126/ synteny system with application to plant genomes. Nucleic Acids Research 39 science.aar7191 (10), 1 May (2011), e68 doi: 10.1093/nar/gkr123 Jiao, W.-B., Accinelli, G. G., Hartwig, B., Kiefer, C., Baker, D., Severing, E., et al. Walker, B. J., Abeel, T., Shea, T., Priest, M., Abouelliel, A., Sakthikumar, S., et al. (2017a). Improving and correcting the contiguity of long-read genome (2014). Pilon: an integrated tool for comprehensive microbial variant detection Frontiers in Plant Science | www.frontiersin.org 9 January 2020 | Volume 10 | Article 1815 Castanera et al. An Improved Melon Reference Genome and genome assembly improvement. PloS One 9, e112963. doi: 10.1371/ loci influencing agronomic traits. Nat. Genet. 51, 1607–1615. doi: 10.1038/ journal.pone.0112963 s41588-019-0522-8 Xiong, Y., and Eickbush, T. H. (1990). Origin and evolution of retroelements based upon their reverse transcriptase sequences. EMBO J. 9, 3353–3362. doi: Conflict of Interest: The authors declare that the research was conducted in the 10.1002/j.1460-2075.1990.tb07536.x absence of any commercial or financial relationships that could be construed as a Ye, K., Schulz, M. H., Long, Q., Apweiler, R., and Ning, Z. (2009). Pindel: a pattern potential conflict of interest. growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads. Bioinformatics 25, 2865–2871. doi: Copyright © 2020 Castanera, Ruggieri, Pujol, Garcia-Mas and Casacuberta. This is 10.1093/bioinformatics/btp394 an open-access article distributed under the terms of the Creative Commons Attri- Zhang, L., Hu, J., Han, X., Li, J., Gao, Y., Richards, C. M., et al. (2019). A high- bution License (CC BY). The use, distribution or reproduction in other forums is quality apple genome assembly reveals the association of a retrotransposon and permitted, provided the original author(s) and the copyright owner(s) are credited and red fruit colour. Nat. Commun. 10, 1494. doi: 10.1038/s41467-019-09518-x that the original publication in this journal is cited, in accordance with accepted Zhao, G., Lian, Q., Zhang, Z., Fu, Q., He, Y., Ma, S., et al. (2019). A comprehensive academic practice. No use, distribution or reproduction is permitted which does not genome variation map of melon identifies multiple domestication events and comply with these terms. Frontiers in Plant Science | www.frontiersin.org 10 January 2020 | Volume 10 | Article 1815