The phased telomere-to-telomere reference genome of Musa acuminata, a main contributor to banana cultivars

Musa acuminata is a main wild contributor to banana cultivars. Here, we reported a haplotype-resolved and telomere-to-telomere reference genome of M. acuminata by incorporating PacBio HiFi reads, Nanopore ultra-long reads, and Hi-C data. The genome size of the two haploid assemblies was estimated to be 469.83 Mb and 470.21 Mb, respectively. Multiple assessments confirmed the contiguity (contig N50: 16.53 Mb and 18.58 Mb; LAI: 20.18 and 19.48), completeness (BUSCOs: 98.57% and 98.57%), and correctness (QV: 45.97 and 46.12) of the genome. The repetitive sequences accounted for about half of the genome size. In total, 40,889 and 38,269 protein-coding genes were annotated in the two haploid assemblies, respectively, of which 9.56% and 3.37% were newly predicted. Genome comparison identified a large reciprocal translocation involving 3 Mb and 10 Mb from chromosomes 01 and 04 within M. acuminata. This reference genome of M. acuminata provides a valuable resource for further understanding of subgenome evolution of Musa species, and precise genetic improvement of banana.


Background & Summary
The wild relatives of domesticated crops, i.e. crop wild relatives (CWRs), generally possess genetic diversity helpful in developing more productive and resilient crop varieties, thereby providing a wide practical gene pool for genetic improvement of crops 1 .In order to address the challenges and threats posed by emerging diseases and climate change, CWRs appear to be a source for solutions to manage both biotic and abiotic stresses 2,3 .At present, combining huge sequence information and precise gene-editing tools provides a route to transform CWRs into ideal crops 2 .Therefore, a high-quality reference genome of CWR germplasm is an important prerequisite for efficiently introducing potential useful genes into breeding programmes.Thanks to the advances in sequencing technologies and analytical tools, many high-quality reference genomes for crops as well as their important wild relatives have been generated.These genetic resources will thus facilitate the identification of structural variants and incorporation of the variants from CWRs into crop gene pools.
Banana domestication started at least 7000 years ago in Southeast Asia 4 .Hybridization between various species and subspecies of the Musa genus led to the development of modern bananas with high production 5 .To date, most banana cultivars were derived from Musa acuminata (A genome), a complex of subspecies geographically segregated in distinct Southeast Asian continental regions and islands 6 .Four particular M. acuminata subspecies have been raised as the main contributors of edible banana cultivars, which are banksii, burmannica, malaccensis, and zebrina 4 .Several large structural variants in these subspecies were identified and suggested to be associated with the domestication of banana [7][8][9][10][11] .Genome research first started in the subspecies malaccensis.The first draft genome of M. acuminata ssp.malaccensis was assembled by incoporating Sanger and Roche/454 reads, with sequence errors corrected by Illumina data 12 .This assembly was anchored along the Musa linkage groups of the genetic map built with SSR and DArT markers.The double-haploid genotype (DH-Pahang) was used in this study for reducing genome complexity and facilitating assembly process.Recently the telomere-to-telomere (T2T) reference genome of DH-Pahang has been constructed using Nanopore data and polished with Nanopore and Illumina reads, with continuity improved significantly 13 .Although DH genotype could miss some important genetic information, these genome resources have significantly facilitated the studies of banana domestication and genome evolution.With advances in the sequencing technologies and biosoftwares, heterozygosity would not be the consistent hurdle.Currently, more and more haplotype-resolved and T2T genomes have been published, such as lychee 14 and apple 15 , providing unprecedented insights into subgenome evolution and domesticated history.
In this study, we assembled a haplotype-resolved and telomere-to-telomere reference genome of M. acuminata ssp.malaccensis by incorporating PacBio HiFi reads, Nanopore ultra-long reads, and high throughput chromatin conformation capture (Hi-C) paired reads.An unphased reference genome was first constructed and used for guiding haplotype-resolved scaffolding (Fig. 1).Multiple assessment methods were applied to evaluate the quality of the haplotype-resolved assembly.A comprehensive genome comparison between this assembly and the previous reference of the DH genotype identified a large reciprocal translocation involving 3 Mb and 10 Mb from chromosomes 01 and 04.Furthermore, the 3-Mb segment (34,734,628 to 37,810,715 bp in chromosome 04) was suggested to be associated with flower development pathway, such as anther/stamen development.The haplotype-resolved genome of M. acuminata will help to obtain a better understanding of potential structural variants, allele specific expression and subgenome evolution of Musa species, and serve as reliable reference for banana breeding programmes.

Methods
Sample collection and sequencing.The M. acuminata sample used for DNA and RNA extraction was obtained from South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, China.Tissues were immediately frozen in liquid nitrogen and preserved at −80 °C for DNA/RNA extraction.The CTAB method was used to extract high quality genomic DNA from leaf tissue samples.

HiFi
Fig. 1 The workflow of generating haplotype-resolved and telomere-to-telomere reference genome for M. acuminata.The unphased reference genome was constructed following the workflow on the left yellow panel.Then a haplotype-resolved and telomere-to-telomere reference genome was produced according to the scheme on the right blue panel.Green boxes represent raw sequencing data; white boxes represent the tools used in this pipeline; pink boxes represent intermediate data; blue boxes represent post-analysis process.
A standard SMRTbell library was constructed using SMRTbell Express Template Prep Kit 2.0 according to the manufacturer's recommendations (Pacific Biosciences, CA, USA) and sequenced on a PacBio Sequel II platform.This yielded 32.39 Gb HiFi data, covering ~65 × coverage of the haploid genome size.The N50 length of the HiFi reads was 17.32 kb.A nanopore library was constructed with the Oxford Nanopore SQK-LSK109 kit following the manufacturers' instructions and sequenced on a PromethION platform.Totally 20.80 Gb ONT data were obtained, covering ~42 × coverage of the haploid genome size.The N50 length was 86.86 kb.A Hi-C library was constructed based on cross-linked genomic DNA and sequenced on an Illumina NovaSeq platform (Illumina, San Diego, CA, USA).In total, 134 Gb Hi-C data were obtained, covering ~268 × coverage of the haploid genome size.The 15.58 Gb NGS data were obtained using the Illumina NovaSeq platform, covering ~31 × coverage of the haploid genome size (Table 1).
Additionally, total RNA was extracted from four tissues, including root, leaf, flower, and fruit, using the NEBNext ® Ultra ™ II Directional RNA Library Prep Kit for Illumina ® (New England Biolabs, MA, USA).
Paired-end 150-bp reads were also generated by the Illumina NovaSeq platform.These yielded a total of 26.90 Gb raw RNAseq data (Table 1).All sequencing were carried out at Anhui Double Helix Gene Technology Co., Ltd.(Anhui, China).
Genome size and heterozygosity estimation.CCS software (https://github.com/PacificBiosciences/ccs)with default parameters was used to generate the consensus reads (HiFi reads).Based on the obtained high-accurate HiFi reads, the K-mer distribution was analysed with jellyfish 16 with jellyfish count -C -m 21 -s 100000000 and jellyfish histo -h 1000000.The results were subsequently imported to GenomeScope v2.0 17 with K-mer length = 21 and Ploidy = 2.The genome size of M. acuminata was estimated to be 450.43 Mb with the 21 K-mer, about 14% shorter than DH-Pahang genome size (523.00Mb) estimated by flow cytometry 12 .
The heterozygosity rate was estimated to be 0.59% (Fig. 2).De novo haplotype-resolved genome assembly.Fastp v0.23.2 18 was performed to filter Hi-C reads with default parameters.Subsequently, hifiasm v0.16.1-r375 19 was carried out to generate the primary unphased draft genome based on HiFi and Hi-C reads.This generated a 491.54Mb draft genome with an N50 of 26.62 Mb, and only 20 contigs consisted of 90% length of the genome (Table 2).Then, ragtag v2.1.0 20with default parameters was first used to sort, orientate, and cluster the primary contigs guided by the T2T version of M. acuminata ssp.malaccensis DH-Pahang genome 13 (Hereafter MAv4).Meanwhile, the primary contigs were anchored into 11 pseudo-chromosomes using Juicer v1.6 21 and 3D-DNA v180922 22 in turn.Then, based on the assembly file obtained from ragtag and the hic file from Juicer and 3D-DNA, Juicebox v2.20.00 23 was introduced for visualizing Hi-C data and manual correction in order to obtain a high-quality reference genome.Finally, there were only 17 gaps in the high-quality reference genome.For gap filling, ONT assembly was constructed by NextDenovo (https://github.com/Nextomics/NextDenovo)with read-cutoff = 1k and genome_size = 500 M. Then this draft ONT assembly was polished by Nextpolish 24 based on the HiFi reads and the Illumina reads with default parameters.Subsequently minimap2 v2.24-r1122 25 with default parameters was used to map the polished ONT assembly to the primary reference genome.We examined the breakpoint with the Integrative Genomics Viewer (IGV) tool 26 and manually filled the gaps based on the alignment results.After using ONT assembly to fill all remaining gaps, a high-quality reference genome named MA was generated.The genome size of this unphased assembly   S1).
To obtain a haplotype-resolved genome, a similar pipeline was applied (Fig. 1).Two primary haploid assemblies were first generated by hifiasm.Further genome assembly statistics were performed with QUAST 27 with default parameters.Accumulative lengths of the two haploid assemblies were 500.78Mb and 484.36 Mb with N50 of 16.53 Mb and 18.58 Mb, respectively (Table 2).After Hi-C scaffolding processes, 469.83 Mb and 470.21Mb were anchored to 11 chromosomes respectively, with an anchored rate of 93.82% and 97.08% (Table 3).The genome sizes of the two haploid assemblies were slightly longer than that of MAv4 (468.82Mb) 13 , and represented approximately 90% of DH-Pahang genome size (523.00Mb) estimated by flow cytometry 12 .All 66 gaps in the two haploid assemblies were filled.Finally, the haplotype-resolved and telomere-to-telomere reference genome for M. acuminata was obtained; and the two haploid assemblies were named MAH1 and MAH2.The circos 28 software was introduced to draw the genome features shown in Fig. 3.The Hi-C heatmap confirmed this assembly as a complete and reliable haplotype-resolved reference genome (Fig. 4).

Genome quality assessment.
Multiple methods were combined to evaluate the quality of genome assembly.First, the HiFi, Illumina, and RNAseq reads were aligned to the phased genome using minimap2 v2.24-r1122, BWA v0.7.17-r1188 29 , and HiSAT2 v2.2.1 30 with default parameters, respectively.BamTools v2.5.1 31 was used to calculate the read mapping rates.The results showed a HiFi coverage rate of 99.86% and 99.87% on MAH1 and MAH2 assemblies, respectively.The mapping rate of Illumina reads reached up to 99.98% in both haploid assemblies.The mapping rate of RNAseq reads ranged from 92.44% to 97.34% (Table 4).Second, the LTR Assembly Index (LAI) calculated from LTR_retriever v2.9.0 32 was used to assess the genome assembly quality.The LAI of MAH1 and MAH2 reached up to 20.18 and 19.48, respectively, indicating that our phased assembly reached the standard of a golden reference.Third, the completeness of the haplotype-resolved genome was evaluated by BUSCO v5.4.3 33 against the 'embryophyta_odb 10' database.In total, 98.57% (1,591 of 1,614) of the complete BUSCO genes were identified (Table 5).Finally, the consensus quality value (QV) of the genome was assessed by Merqury v1.3 34 with meryl k = 19 count, showing 45.97 and 46.12 of QV (Genome accuracy >99.99%) for MAH1 and MAH2, respectively (Table 6, Supplementary Figure S2).S1).These results were comparable with the findings in the previous T2T DH genome version (Repeat elements: 52.62%; LTR: 34.85%) 13 .Standard MAKER3 v3.01.03 38 pipeline was used to annotate genes.All high-confidence protein sequences in swiss-prot 39 database were imported for homology prediction.Transcripts from the 4 tissues, including root, leaf, flower and fruit, were used for gene prediction.Then AUGUSTUS v3. 3  train the ab-initio gene models.Finally, the MAKER3 pipeline was run again to obtain high-quality gene annotations.Functional characterization of the predicted coding genes was performed using eggNOG-mapper v2 40 based on the eggNOG v5.0 database 41    (74.72%) genes were functionally annotated (Table 7).Besides, BUSCO scores of protein-coding genes in MAH1 and MAH2 were up to 89.41% and 90.27% (Table 8).
Identification of telomeres and centromeres.TIDK v0.2.1 (https://github.com/tolkit/telomeric-identifier) was used to find telomeres.In total 36 telomeres were found (Table 9).Plant centromeric regions are generally characterized by the presence of short tandem repeats that are highly enriched in these regions 42 , accompanied by a collapse in the density of LTR elements such as Copia.By identifying these distinctive features, centromeric regions can be located.We predicted centromeric regions according to the workflow in Shi et al. 43 , which employed the above approach.Using Tandem Repeats Finder v4.09 44 with the parameters: trf genomes.fa 2 7 7 80 10 50 500 -f -d -m, we screened 183 bp, 148 bp, 124 bp, 125 bp, and 191 bp tandem repeat units as candidates based on sorted results and IGV results (Supplementary Table S2, Supplementary Figure S3).The centromeric regions were defined according to the density of 183 bp tandem repeat unit, which was the highest enriched centromeric repeat unit.Finally, all centromeric regions have been captured successfully (Table 10, Supplementary Figure S3).
Characterization of a reciprocal translocation in Musa acuminata.Nucmer v4.0.0rc1 45 was used to obtain the syntenic relationship between MAH1 and MAH2 with default parameters.Then the delta-filter was launched with parameters '-i 90 -l 15000' .In the same way, our haplotype-resolved assembly was aligned against MAv4 using nucmer.Mummerplot command was used to generate the dot plots (Supplementary Figure S4).Syri v1.6.3 46 with default parameters was used for identifying structural variants between MAH1 and MAH2 (Fig. 5).Overall, 47 translocations with a cumulative size of 2.70 Mb (~0.57%), 23 inversions with a cumulative size of 11.30 Mb (~2.40%), and 53 duplications with a cumulative size of 1.33 Mb (~0.28%) were defined.These structural variants were generally heterozygous, representing more complete genetic information compared with the double-haploid MAv4 genome.MCscan tools 47 were used to search for the syntenic relationships between the two haploid assemblies and MAv4 at the gene level.Briefly, 'jcvi.compara.catalog'module with '--cscore = 0.99' and 'jcvi.compara.synteny'module with '--minspan = 30' were used to build the syntenic regions; then, syntenic relationships were visualized by 'jcvi.graphics.karyotype'module.Besides, potential structural variants and heterozygous regions were shown in Supplementary Figure S5.A reciprocal translocation involving 3 Mb and 10 Mb from chromosome 01 and 04 was identified (Fig. 6a).These reciprocal translocation gene blocks were located in the translocated regions identified in whole genome alignment results (Supplementary Figure S4C,D).The 10-Mb segment from 261,650 to 10,745,936 bp in chromosome 01 of MAH1 was linked to 44,882,868 to 34,419,170 bp in chromosome 04 of MAv4 (Supplementary Figure S5).The 3-Mb segment from 34,734,628 to 37,810,715 bp in chromosome 04 of MAH1 was linked to 122,362 to 3,101,126 bp in chromosome 01 of MAv4.The reciprocal translocation between MAH2 and MAv4 was located in the similar genomic regions.The huge difference in chromosome length in chromosome 01 and chromosome 04 between MAH1/2 and MAv4 was also derived from this reciprocal translocation, while other chromosome lengths and genomic total lengths were comparable (Table 3).
We further performed GO enrichments based on the extracted genes located in the translocated regions using TBtools v1.108 48.The genes in the 10-Mb segment of MAH1 were not significantly enriched in any biological process, while those on the 3-Mb segment were enriched in several pathways associated with flower development (Supplementary Table S3), including anther development (GO:0048653), stamen development (GO:0048443), regulation of flower development (GO:0009909), and floral whorl development (GO:0048438).For further validation, we used nucleotide BLAST tools in National Center for Biotechnology Information Fig. 5 The sequence collinearity and structural variants between MAH1 and MAH2.
(NCBI) with default parameters to align identified genes to non-redundant database, and checked gene functions manually.

Data Records
All raw sequencing reads have been deposited in the National Center for Biotechnology Information (NCBI) under BioProject ID PRJNA962682, and the National Genomics Data Center (NGDC) under BioProject ID PRJCA018611.The PacBio HiFi, Nanopore, Hi-C, Illumina sequencing reads have been deposited in the NCBI Sequence Read Archive database with accession group numbers SRP435127 49 .Genome assembly is available from GenBank with accession number GCA_030219345.1 50 .The genome annotation files have been submitted to the online open access repository Figshare database 51 , including a high-quality reference genome that we constructed for guidance.

Technical Validation
Manual correction for chromosome scaffolding.For constructing a high-quality reference genome, we used Juicebox to manually correct the reference genome based on Hi-C alignments (Supplementary Figure S1).Finally, 471.04 Mb (95.83%) contigs were anchored to 11 pseudo-chromosomes.Then, we started to orient, sort and group our haplotype-resolved genome based on this high-quality reference genome.We also used Juicebox to manually correct the haplotype-resolved genome based on Hi-C alignments.In total, 469.83 Mb (93.82%) and 470.21Mb (97.08%) contigs were anchored to 11 chromosome pairs, respectively.We further examined the Hi-C alignments in chromosome 01 and chromosome 04 in Juicebox (Fig. 6b), and confirmed the accurate assemblies of chromosome 01 and chromosome 04.Besides, chromosome 01 consists of only one contig (Table 3), further confirming its high continuity.
Targeted pCr confirmed the reciprocal translocation between Chr01 and Chr04.Based on the genomic syntenic analysis between our assembly and MAv4, we identified a large reciprocal translocation from chromosomes 01 and 04, corresponding to the translocation found in a previous study 9 .In that study, three pairs of primers were designed to amplify the breakpoints located along the reference and hypothesized chromosome structures, thereby showing the presence of chromosomes 01, 04, and 1T4 resulting from the translocation.Here we used the same primer pairs to perform targeted PCR to validate the chromosome structures found in our sample (Fig. 6c).DNA was extracted from leaf tissue of M. acuminata ssp.malaccensis.PCR was performed in 50-μL volumes containing 2.5 ng of gDNA, 1 μL of specific primers, 32 μL of distilled, deionized water, and 0.5 μL of TaKaRa LA Taq ® (Vazyme) using an Eastwin Life Science EDC810 PCR amplification system.The reaction conditions for thermal cycling were 94 °C for 5 min, followed by 35 cycles of 94 °C for 45 s, 56 °C for 45 s, and 72 °C for 60 s.Thereafter, PCR products were visualized by 2% agarose gel-electrophoresis with a 100 bp DNA ladder.Only the breakpoint of chromosome 1T4 was amplified in our studied sample, suggesting that the reciprocal translocation involving 3 and 10 Mb segments from chromosomes 01 and 04 existed in both haploid genomes of the M. acuminata sample (Fig. 6c).This finding was consistent with our whole genome alignment results in sequence and gene levels.

Fig. 3
Fig. 3 The Overview of M. acuminata genome assembly and features.The tracks represent the following elements (from outer to inner): (a) Karyotypes of the 22 chromosome sequences, (b) TRF-183bp centromeric repeat density, (c) Copia density, (d) Transposable element (TE) density, (e) GC contents, (f) Gene density.The innermost is syntenic relationships.

Fig. 6
Fig. 6 The large reciprocal translocation between chromosome 01 and chromosome 04.(a) Genomic synteny between DH-Pahang (MAv4) and wild M. acuminata.(b) Hi-C heatmap of chromosome 01-04.(c) The possible chromosome structures in wild M. acuminata.The targeted PCR result suggests that homozygous reciprocal translocations (only 1T4) exist in our sample, consistent with our whole genome alignment results in sequence and gene levels.

Table 1 .
Summary of sequencing data of Musa acuminata ssp.malaccensis for haplotype-resolved and telomere-to-telomere assembly and genome annotation.

Table 2 .
Summary of genome assembly of Musa acuminata ssp.malaccensis genome.Note: MA represents the primary contig sets, while MAH1 and MAH2 represent contigs in haplotype1 and contigs in haplotype2.

Table 3 .
The lengths of the pseudo-chromosomes of Musa acuminata ssp.malaccensis genomes. is 471.04Mb with an anchored rate of 95.83%.The Hi-C heatmap confirmed the contiguity of the assembly (Supplementary Figure

Table 4 .
Assessment of genome quality based on mapping with RNAseq reads.

Table 5 .
BUSCO homology repeats in RepeatMasker.In total we identified 235.46 Mb (50.11%) and 234.61Mb (49.90%) repetitive sequences in MAH1 and MAH2, respectively.Among these, long terminal repeats (LTR) that accounted for 36.61% in MAH1 and 34.19% in MAH2 were the most abundant repeat elements (Supplementary Table

Table 6 .
.2 and SNAP v20131129 were used to The consensus quality values of MAH1 and MAH2.

Table 7 .
Statistics of protein-coding genes in MAH1 and MAH2.

Table 9 .
Summary of telomere information of Musa acuminata ssp.malaccensis genome.

Table 10 .
Summary of centromere information of Musa acuminata ssp.malaccensis genome.