In order to obtain the T2T A. eriantha genome, we introduced multiple sequencing technologies to generate high-quality sequences (Supplementary Table 1). PacBio sequencing yielded 29.7 Gb of HiFi clean reads (~ 49 X genome coverages) with an N50 of 17.1 Kb; Oxford Nanopore Technologies yielded 20.6 Gb (~ 34X genome coverages) with an N50 of 100.4 Kb; A total of 105.5 Gb of clean reads (~ 175X genome coverages) were generated from Hi-C libraries; The Illumina short-reads data (~ 8X genome coverages) from “Midao31” were used to polish ONT long-reads.

First, the HiFi reads were used for primary assembly of A. eriantha by hifiasm, generating two preliminary assemblies with a size of 633 Mb and 622 Mb. The two preliminary assemblies include 237 contigs with the N50 length of 21 Mb and 185 contigs with the N50 length of 18 Mb respectively, which are greatly larger than the previous two versions of A. eriantha (Table 1). Then the contigs of two haplotypes were anchored to 29 scaffolds by juicer (Durand et al. 2016) and 3D de novo assembly (3D-DNA) pipeline (Dudchenko et al. 2017) using Hi-C reads. It was worth noting that the 29 scaffolds still contain a lot of gaps and incorrect joins because the HiFi contigs were optionally split in regions lacking Hi-C coverage. Then we re-mapped the HiFi contigs of two haplotypes against the corresponding 29 Hi-C scaffolds for de novo genome assembly with a reference-guided strategy by using a custom Perl script (https://github.com/aaranyue/CTGA). We consequently obtained two preliminary chromosome-level genomes that contain 8 gaps and 17 gaps, respectively. Meanwhile, the polished ONT reads were used to fill gaps using the TGS-GapCloser (Xu et al., 2019). Finally, we get two final haplotype A. eriantha genomes, termed MDHAPA and MDHAPB, with a size of 619.3 Mb and 611.7 Mb, respectively, including 27 and 28 gap-close chromosomes (Table 1). The remaining three gaps were distributed on chr13, chr21 of MDHAPA, and chr21 of MDHAPB (Supplementary Table 2). Our two haplotype assemblies, with greatly improved contiguity and completeness, showed high synteny with the four published genomes of “White”, “wild”, “Hongyangv4.0” and “Red5” (Supplementary Fig. 2).

Within the two haplotype assemblies, 253.46 Mb (40.92%) and 251.48 Mb (41.11%) of repetitive sequences were found (Table 1), similar to those in the genomes of “White” (43.32%) and “wild” (41.29%). Among them, 25.45% and 24.67% were long terminal repeats (LTRs), and 15.24% and 12.45% were terminal inverted repeats (TIRs) in MDHAPA and MDHAPB, respectively (Supplementary Table 3). We predicted a total of 46,008 and 47,184 high-confidence protein-coding genes from MDHAPA and MDHAPB, respectively, which were close to the 45,809 genes predicted in the HY4P genome and the 45,434 predicted in the HY4A genome (Yue et al. 2022) and showing 93.8% and 93.2% BUSCO completeness (Table 1). The total lengths of the MDHAPA and MDHAPB genes were 249.73 Mb and 252.05 Mb, respectively. Of these protein-coding genes, 39,836 (86.58%) and 40,328 (85.47%) could be functionally annotated to a comprehensive database of eggNOG. Compared to the “White” genome, more than 4000 genes were newly predicted in the two assemblies. And a large number of genomic variations between MDHAPs and White were discovered to be associated with these new predicted genes (Supplementary Table 4). KEGG analysis revealed that these genes functionally enrich in multiple pathways, including RNA degradation, plant-pathogen interaction, pentose glycosylation, and so on (Supplementary Fig. 3a and 4a). GO analysis showed that the newly predicted genes were mainly enriched in the biological process and molecular function categories (Supplementary Fig. 3b and Fig. 4b). Among them, a lot of genes were related to cold acclimation and GTPase activity. So, these newly predicted genes may involve important biological functions in kiwifruit.

Multiple evaluative approaches were used to confirm the high quality of haplotype assemblies: (i) The assembly spectrum copy number plots plotted by the KAT program (Mapleson et al. 2017) clearly demonstrated that the phasing of the assembled haplotypes is correct (Supplementary Fig. 5), (ii) The 97.25% Illumina reads, the 99.49% HiFi reads, and the 99.99% were successfully aligned to the two assemblies (Supplementary Table 5), suggesting that the A. eriantha genome was sufficiently covered by the assembly; (iii) Hi-C interaction matrices displayed a diagonal pattern for the intra-chromosomal interactions in all pseudochromosomes, indicating the right ordering and orientation (Fig. 1b and 1d); (iv) A total of 1597 (99.0%) and 1601 (99.2%) complete BUSCO genes were identified in the two haplotype genomes, respectively, which was superior to the genome of the “White” and “wild” (Supplementary Fig. 6); (v) The high LAI scores (21.52 and 20.45) were observed in two haplotype genomes (Table 1). (vi) The consensus quality value (QV) of MDHAPs was 59.60 and 50.93, indicating the high accuracy of our assemblies (Table 1).

Syntenic analysis revealed highly consistent sequence orders between the two haplotypes (Fig. 2a). Genomic collinearity analysis between MDHAPA and MDHAPB revealed that 45,295 transcripts matched 540 syntenic blocks (Fig. 2b). We also identified substantial variations between the two haplotypes by using SyRI, including 2,574,323 SNPs, 457,407 insertions/deletions, 289 inversions, and 167 translocations (Supplementary Table 6). These variations spanned approximately 35.7 Mb, representing 5.6% and 5.7% of the assembled MDHAPA and MDHAPB genomes, respectively.

Correlated or differential expression of alleles could have profound effects on growth and evolution. Using MCscanX, we separated 35,395 pairs of alleles in the haplotype-resolved genome (Fig. 2c). These alleles are relatively evenly distributed across the 29 chromosomes of kiwifruit. Most allelic genes maintained high levels of coding sequence similarity (mean = 98.13%) (Fig. 2d). To evaluate the natural selection pattern of alleles, we calculated the K_a/K_s value between allelic gene pairs. The results revealed that the majority of allelic genes had clearly experienced purifying selection (K_a/K_s < 1), with only a small number of allelic pairs showing possible positive selection (K_a/K_s > 1) (Fig. 2e).

Then we investigated the allele-specific expression (ASE) across different tissues (fruit at different development stages, leaf and stem). A total of 2904 allelic genes showed different expression patterns (log2 fold change > 1, p value < 0.05) (Fig. 2f), while others did not show a more than two-fold difference in expression, indicating that most alleles in the A. eriantha genome were coordinately expressed. These genes showed functional enrichment in biological processes such as regulation of lipids, activation of protein kinase activity, and cellular components including polysomal and vesicle tethering complexes based on the gene ontology (GO) annotation (Supplementary Fig. 7a). The result of KEGG pathway annotation revealed that the ASEs had functional enrichment in multiple biological processes, including endocytosis, cellular senescence, and circadian rhythm (Supplementary Fig. 7b), suggesting that a potential mechanism to overcome deleterious mutations occurred in important genes related to basic biological functions.

Among the 2904 ASEs, 786 genes showed expression biased toward one allele across six samples, which were defined as consistent ASE genes (Supplementary Fig. 8a). Genes with expression biased toward one parental allele in some samples but shifted to another allele in other samples were defined as inconsistent ASE genes, which indicate an overdominance effect. A total of 103 inconsistent ASE genes were found in the A. eriantha “Midao31” genome (Supplementary Fig. 8b), which was less than the number of consistent ASE genes. Then we investigated the functional impact of the allelic variations and found that 29.97% of SNPs and 30.88% of InDels caused changes in the upstream sequence (+ 2 kb) (Supplementary Table 7), indicating that the promoter sequence variation may be an important cause of allelic differential expression.

Telomeres are nucleoprotein structures at the ends of chromosomes and function to maintain genome stability and consist of a tandem repeat of TG-rich microsatellite sequence (Turner et al. 2019). The telomeric repeat sequenece is TTTAGGG in most plants (Fajkus et al. 2005). Using the seven-base telomeric repeat (CCCATTT at the 5' end and TTTAGGG at the 3' end) as a query, we identified 54 telomeres, resulting in 25 T2T pseudomolecules both in MDHAPA and MDHAPB assemblies (Supplementary Table 8 and Supplementary Table 9). The total counts of telomere repeats ranged in 176 ~ 3140 and 111 ~ 2760, with a mean value of 1318 and 1063.

The centromere is an important part of chromosomes and plays a crucial role in the proper segregation of chromosomes. Unfortunately, information about the centromeres of kiwifruit is limited. One of the main reasons is that the centromere region contains highly repetitive sequences, which impede assembly from short DNA sequencing reads (Nurk et al. 2022). Nevertheless, the development of sequencing technology allow us to have the opportunity to assemble centromere regions. To identify the location and sequence features, we used the Tandem Repeats Finder (TRF) tool to search tandem repeats in our assemblies, and only the repeat monomers with lengths ranging from 100 to 200 bp were retained. And then CD-HIT (Fu et al. 2012) was used for clustering these monomers to reduce sequence redundancy and improve the precision of centromere localization based on sequence similarity search, the continuous and high-frequency regions were thought to be approximate centromeric sequences. Finally, we determined the location of the centromeres of all chromosomes in two haplotype assemblies. The result showed that the centromere boundaries of the two haplotype genomes had similar positions on the chromosome, and the length of the centromere region ranged from 217,369 bp to 1,893,971 bp in MDHAPA and from 112,182 bp to 1,168,845 bp in MDHAPB (Supplementary Table 10). Aside from that, there are 147 and 151 new genes predicted in the centromere region of MDHAPA and MDHAPB, respectively (Supplementary Table 11 and 12). To verify the accuracy of the centromere region, we analyzed the gene density, repeat distribution, and sequence similarity on the chromosome (Fig. 3a). The distribution of repeats revealed that the class I retrotransposons are more common in centromeres, while the class II retrotransposons were more evenly distributed across the genome, which was similar to other species such as Brassica (Perumal et al. 2020) and A. chinensis (Yue et al. 2022). In addition, the centromere region has low gene density and low similarity compared with other regions on the chromosome (Fig. 3a). Finally, the Hi-C heatmap also showed that the location of centromere region was correct(Fig. 3b).

Based on the analysis of tandem repeats in the centromere region, we found that the repetition types in the centromere region of kiwifruit are quite complex, containing one major repeating monomer and several minor repeating monomers. The length of the major repeating monomer ranged from 102 to 192 bp (Supplementary Table 13 and 14), and most of the major repeating monomer’s length was 153 bp (named Ae-CEN153) (Fig. 3c), which was consistent with A. chinensis. Alignment of the Ach-CEN153 and Ae-CEN153 sequences from A. chinensis and A. eriantha. Alignment of the Ach-CEN153 and Ae-CEN153 sequences from A. chinensis and A. eriantha revealed that they retain strong similarities (Supplementary Fig. 9). These results indicated that the centromere repeat unit sequences were relatively conserved in the two species. The phylogenetic analysis of the representative monomers of rice, Arabidopsis, A. chinensis and all monomers of A. eriantha revealed that all monomers could divided into four groups, the kiwifruit monomers were separated from Arabidopsis and rice (Fig. 3d), implying that the kiwifruit centromere tandem repeats have little homology with rice and Arabidopsis.

“Midao 31”, a hybrid between A. eriantha “White” (female) and A. eriantha “MHX-1” (male), was used in this study. High-quality genomic DNA was extracted from fresh young leaf tissue from “Midao 31”, growing in Hefei, Anhui Province, China, and separately packaged for PacBio HiFi, ONT ultra-long, and Hi-C sequencing. Tissue materials from leaves, stems, and fruits were used for RNA-seq and genomic annotation.

The kiwifruit genomic DNA (Allen et al. 2006) was prepared using a modified CTAB method and evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and Qubit fluorometer instrument (Thermo Fisher Scientific, MA, USA). The library for ONT sequencing was constructed using the 1D ligation sequencing kit (SQK-LSK108, ONT, UK) and sequenced on the ONT PromethION platform. The Hi-C libraries were then constructed as a standard procedure, including chromatin extraction and digestion, DNA ligation, and purification. For PacBio HiFi sequencing, a standard SMRTbell library was prepared with 50 μg of gDNA by using the SMRTbell Express Template Prep Kit 2.0, according to the manufacturer’s instructions. SMRTbell libraries were then sequenced on a PacBio Sequel II system (Pacific Biosciences, CA, USA).

The HiFi reads were de novo assembled into two haplotype contigs using the hifiasm v0.16.1 software with default parameters (Cheng et al. 2021). And then the two haplotype contigs were corrected, grouped, sorted, and anchored to pseudochromosomes with Hi-C reads using juicer (Durand et al. 2016) and 3D-DNA software (Dudchenko et al. 2017). Then, a custom Perl script was used to remap the HiFi contigs against the Hi-C pseudochromosomes using a reference-guided strategy (https://github.com/aaranyue/CTGA). The ONT reads were polished by Pilon (Walker et al., 2014) and then used to fill gaps in the ref-guided pseudochromosomes by the TGS-GapCloser (Xu et al., 2019). To evaluate the reliability of our assembly, the HiFi reads and the Illumina reads were remapped to the two assembled haplotypes genomes using minimap2 (Li et al., 2018). The completeness of two assemblies was evaluated by mapping the Benchmarking Universal Single-Copy Orthologs (BUSCO) to the genomes using BUSCO v3.0.2 (Manni et al. 2021) with the Embryophyta odb10 dataset. The continuity was assessed using the LTR Assembly Index (LAI) (Ou et al., 2018). The consensus quality (QV) value was assessed using Mequery (v1.3) (Rhie et al. 2020).

Gene structure annotation was performed by the braker software (Hoff et al. 2019) using a combination of de novo prediction and transcript evidence from ten RNA-seq datasets (Supplementary Table 15). Only genes that met the criteria of having start and stop codons and being longer than 100 nucleotides were reserved. Gene function was annotated by eggmapper against a series of protein sequence databases. GO and KEGG enrichment analyses were performed using the R package clusterProfiler (Yu, 2012). Transposable elements (TEs) in the A. eriantha genome were identified using the Extensive de novo TE Annotator (EDTA) with default parameters (Ou et al. 2019). The tandem repeats (TRs) were identified by the TRF software (Benson 1999) with parameters (2 7 7 80 10 50 500 -f -d -m).

The genomic sequences of “White” and “wild” were downloaded from the KGD database (http://kiwifruitgenome.org/) (Yue et al. 2020) and the NCBI database. The MUMmer software (Marcais et al., 2018) was used for the alignment comparison analysis between MDHAPA, MDHAPB, and the “White” and “wild” genomes with parameters (-maxmatch -c 500 -b 200 -l 100). The alignment was filtered using the delta-filter implemented in Mummer with the parameters (-m -i 90 -l 2000). Then, the filtered result was visualized by mummerplot. The collinearity, structural variations, and sequence differences analysis was performed between MDHAPA and MDHAPB using the Synteny and Rearrangement Identifier (SyRI) (Goel et al. 2019).

MCScanX (Wang et al., 2012) was used to identify synteny blocks between a pair of allelic chromosomes, paired genes within each synteny block with high similarity were considered alleles A and B. The K_a/K_s value was calculated using KaKs_Calculator V2.0 (Wang et al., 2010). The RNA-seq data from 6 samples (fruit DAF 20, fruit DAF 40, fruit DAF 60, fruit DAF 120, leaf, and stem) were aligned to the MDHAPA and MDHAPB using HISAT2 v2.0.0 (Kim et al. 2015). The FPKM values were estimated using featureCounts v1.5.3 (Liao et al. 2014). The differentially expressed genes were identified using DESeq2 (Love MI et al., 2014) and edgeR (Robinson MD et al., 2010). ASE was determined if the log₂(fold change) values of FPKM between two alleles were greater than 1 and the P value < 0.05. The ASE patterns (consistent and inconsistent expression pattern) was identified using the same methods as the tea plant Camellia sinensis genomes project (Zhang et al. 2021).

The telomere was identified by the Telomere Identification toolKit (tidk) using the normalized and unified sequence “AAACCCT” for search. A pipline containing TRF tools, and CD-HIT tools was used for the identification of the centromere region. The TRF tools were used to identify whole-genome tandem repeats and monomers. The CD-HIT was used to cluster these monomers to reduce sequence redundancy. Those continuous and high-frequency regions were regarded as candidate centromere regions. At last, combine the result of gene density, TE number with candidate regions to predict the most likely centromere location.

The representative centromeric monomer sequences of Arabidopsis (Wang et al. 2022), rice (Song et al. 2021), A. chinensis and A. eriantha were then aligned by muscle (Edgar 2004), and phylogenetic analysis was performed by the Maximum Likelihood method with default parameters and then visualized with iTOL (Letunic et al., 2007). The Ach-CEN153 and Ae-CEN153 sequences were aligned and visualized using DNAman V9.0.

Not applicable.

All authors hereby consent to publication of the work.

The authors declare no competing financial interests.