The complete genome assembly of Nicotiana benthamiana reveals genetic and epigenetic landscape of centromeres

doi:10.21203/rs.3.rs-4568898/v1

Download PDF

Article

The complete genome assembly of Nicotiana benthamiana reveals genetic and epigenetic landscape of centromeres

https://doi.org/10.21203/rs.3.rs-4568898/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Nicotiana benthamiana is a model organism widely adopted in plant biology and biotechnology. Its genomic research has lagged since its initial release in 2012. To further improve its usefulness, we generate and phase the complete 2.85 Gb genome assembly of allotetraploid N. benthamiana with all 19 centromeres and 38 telomeres fully resolved. We find that although Solanaceaecentromeres are widely dominated by Ty3/Gypsy retrotransposons, satellite-based centromeres are surprisingly common in N. benthamiana with 11 of 19 centromeres featured by megabase-scale satellite arrays. Interestingly, the satellite-enriched and satellite-free centromeres are extensively invaded by distinct Gypsy retrotransposons where CENH3 protein more preferentially occupies, suggestive of their crucial roles in centromere function. We demonstrate that rDNA is a major origin of centromeric satellites, and mitochondrial DNA could be employed as core component of centromere. Subgenome analysis indicate emergence of satellite arrays likely drives the centromere formation and maturation during genomic shock after polyploidization. Altogether, we propose N. benthamiana centromeres evolves via neocentromere formation, satellite expansion, retrotransposon enrichment, and mtDNA integration.

Biological sciences/Genetics/Genomics/Genome evolution

Biological sciences/Plant sciences/Plant evolution

Centromere is a specialized domain that attaches chromosomes to spindle fibers during cell division, responsible for the sister chromatid cohesion¹. Despite its conserved function, centromeric sequence evolves rapidly, exhibiting remarkable diversity among species^2,3. The centromeric chromatin is defined by the presence of a centromere-specific histone H3 (CENH3) in plants. Centromeres are composed of long arrays of satellite repeat in many eukaryotic organisms, such as human⁴, rice⁵ and Arabidopsis⁶. The centromeric satellite generates higher-order repeats (HORs) through homogenization, resulting in whole centromeres dominated by a single repeat³. Interestingly, not all centromeres contain canonical satellite arrays. The Ty3/Gypsy retrotransposons serving as core of functional centromeres are found in several Poaceae plants including wheat⁷, einkorn⁸, oat⁹ and maize¹⁰, and Solanaceae plants such as potato¹¹, tomato¹², tobacco¹³ and pepper¹⁴. The enriched centromeric retrotransposons contribute to centromere function in a positive manner, contrary to the detrimental role in human⁴ and Arabidopsis⁶ centromeres. The genetic and epigenetic pattern of satellite-based centromeres has been well studied in humans⁴ and Arabidopsis⁶. CENH3 occupies homogenized higher-order repeats and the invaded transposons disrupt the genetic and epigenetic organization^3,6. However, the chromatin landscape of Gypsy-based centromeres is poorly understood.

Nicotiana benthamiana is a model plant widely used in plant biology and biotechnology owing to its fast transient transgene analysis. The N. benthamiana originates from an interspecific hybridization at 5–6 million years ago (Mya) leading to allopolyploid formation, chromosome loss and extensive chromosomal rearrangements¹⁵. As such, an understanding of the adaptation and evolution of centromeres to “genome shock” due to allotetraplodization event requires further study. Complete assembly of centromere is challenging due to the presence of abundant repetitive sequences in and around centromeres. Leveraging advances in assembly algorithm¹⁶ and sequencing technology such as PacBio HiFi and ultra-long ONT (Oxford Nanopore Technology), many complete telomere-to-telomere (T2T) genome assemblies have been reported for human⁴, rice⁵, Arabidopsis⁶ and maize¹⁰. The N. benthamiana genome was first sequenced in 2012; until recently, multiple updates were reported with contig N50 reaching 54.2 Mb, yet the centromeres and telomeres remain incomplete^15,17–20. Generating the complete T2T genome sequence is paramount for fundamental and applied research using N. benthamiana and dissecting its mysterious centromere landscape.

Here we achieve the first complete genome assembly of the allotetraploid model plant N. benthamiana LAB strain using latest long-read sequencing, which represents a major step forward for the polyploid genome assembly. Our analyses unravel the evolutionary history of N. benthamiana and offer insights on the genetic and epigenetic landscape of centromeres, which greatly enrich our knowledge on polyploid centromere diversity and evolution. The complete genome assembly and annotation provide valuable resources that will enhance the application of N. benthamiana as a model organism in plant biology and biotechnology.

Assembly and annotation of the complete genome

We have generated a complete assembly of N. benthamiana using a combination of PacBio HiFi, ONT ultra-long, Hi-C sequencing and Bionano optical mapping. A total of 332.6 Gb PacBio HiFi reads and 136.5 Gb Nanopore ONT reads were generated, representing 116.7× and 47.9× coverage of the N. benthamiana genome, respectively (Supplementary Table 1). The initial genome assembly was performed using hifiasm¹⁶ integrating HiFi and ultralong ONT reads (Supplementary Fig. 1), resulting in a 2,845.07 Mb draft genome containing all 38 telomeres, with a contig N50 length of 144.74 Mb (Supplementary Fig. 2). Then gap filling and NOR (nucleolus organizer region) assembly resolved the remaining three gaps in repetitive regions. One gap related to nuclear plastid DNA segments (NUPTs) was closed by ultralong ONT reads, which was supported by the Bionano molecules (Extended Data Fig. 1a). Two gaps in NOR regions were resolved by adding a length of 3.18 Mb rDNA arrays containing 327 copies of 45S rDNA to Chr06, and a length of 9.62 Mb arrays containing 920 copies of 45S rDNA to Chr19, respectively (Supplementary Fig. 3). After telomere patching, we finally obtained a complete T2T assembly of the N. benthamiana genome, with a total size of 2,849.3 Mb and contig N50 of 146.41 Mb (Fig. 1a,b and Table 1).

Table 1

Statistics for two genome assembly and annotation of *N. benthamiana*
Genomic feature	NB.PCP ¹⁸	NB.T2T
Number of contigs	1,901	19
Total length (bp)	2,926,018,961	2,849,299,002
Contig N50 (bp)	31,273,564	146,412,814
Number of gaps	437	0
Number of telomeres	3	38
Number of centromeres	19	19
Number of gene models	57,583	57,023
GC content (%)	38.12	38.08
Repeat content (%)	78.01	78.79
Assembly BUSCOs (%)	98.80	98.80
Annotation BUSCOs (%)	95.60	99.30
LAI	17.85	17.84
QV	60.23	69.63

The N. benthamiana complete genome (NB.T2T) showed substantial improvements over previous assemblies^15,17–20 (Supplementary Tables 2 and 3), representing the first complete genome of polyploid Nicotiana plants. A total of 50.47 Mb and 144.42 Mb sequences were newly assembled to chromosomes compared to the recently published genome (NB.PCP)¹⁸ and the widely used NB.261 (version 2.6.1, http://solgenomics.net) assemblies, respectively (Fig. 1c). Take the NB.PCP as an example, the majority of the new assembled sequences were in difficult-to-assemble repetitive regions including telomeres and subtelomeres (48.5%), rDNA arrays (30.1%), and segmental duplications (18.5%). We totally annotated 2,247.3 Mb repetitive sequences in the N. benthamiana genome, accounting for 78.87% of the assembly, which contained 63.50% retrotransposons and 5.22% DNA transposons (Supplementary Table 4). LTR retrotransposons occupied 59.60% of the whole genome, including 42.40% Ty3/Gypsy and 16.47% Ty1/Copia elements. In addition, we identified 30.26 Mb microsatellites, 31.02 Mb satellites, 8.71 Mb 5S rDNAs, 13.98 Mb 45S rDNAs, 1.23 Mb NUPTs and 7.44 Mb nuclear mitochondrial DNA segments (NUMTs), which together accounted for 3.25% of the genome.

Genome annotation predicted 57,023 protein-coding genes in NB.T2T, of which 84.23% were functionally annotated (Supplementary Fig. 4). The NB.T2T annotation shared 50,881 genes in 33,926 orthogroups with two recently published versions^18,20, accounting for 89.2% of all annotated protein-coding genes (Supplementary Fig. 5). It is worth noting that the NB.T2T showed higher completeness and consistency in proteome assessment than other versions (Fig. 1d and Supplementary Fig. 5), providing an valuable resource for plant research community.

Extensive validation of the T2T assembly

To verify the primary hybrid assembly, we also assembled the HiFi reads using hifiasm¹⁶ and obtained a HiFi assembly (2,860.66 Mb, contig N50: 92.85 Mb) with 28 gaps (Supplementary Fig. 6). We compared our HiFi assembly with two previously published HiFi assemblies^18,19, and found 21 ‘conserved’ gaps (Extended Data Fig. 1b), among which 17 were at GA/TC low-complexity repeat regions (Extended Data Fig. 2 and Supplementary Table 5). The HiFi chemistry bias is responsible for contig breaks at GA- or GAA-rich regions, highlighting the importance of ONT sequencing in achieving the T2T assembly²¹. The remaining four conserved gaps were located within two 45S rDNA arrays, one segmental duplication (Fig. 2c) and one NUPT insertion (Fig. 2d), which were well closed in the final NB.T2T assembly with high coverage of HiFi or ONT mapped reads spanning these regions.

Then we applied multiple strategies to evaluate the accuracy and completeness of the NB.T2T assembly. Firstly, we manually checked the Hi-C chromatin interaction maps and Bionano optical maps with neither showing any obvious mis-assembly (Fig. 1b and Extended Data Fig. 1a). Then whole-genomic alignment exhibited strong collinearities between each of two recently published genome assemblies of N. benthamiana and NB.T2T (Fig. 1e and Extended Data Fig. 1b), demonstrating the large-scale structural accuracy of NB.T2T. Secondly, mapping Illumina, HiFi and ONT reads back to NB.T2T yielded high mapping rates (> 99.5%), and showed nearly uniform coverage across whole genome (Fig. 2a and Extended Data Fig. 2). The overall base accuracy was 99.95% based on variant calling using HiFi reads (Fig. 2b), confirming the overall accuracy of the base-scale. Thirdly, we determined that NB.T2T had a QV (quality value) of 69.65 and completeness value of 99.14% using Merqury²² (Table 1). LTR assembly index (LAI) of 17.84 reached quality standard for reference genome²³. BUSCO evaluation of the assembly identified 98.81% of the complete genes (Fig. 1d), further demonstrating the high completeness of NB.T2T. Lastly, all telomeres were detected using the telomeric repeat of “AAACCCT” and sub-telomeres using the 181-bp satellite repeat (Extended Data Fig. 3). The genomic synteny and long-read coverages confirmed the complete assembly of telomeric and subtelomeric regions. Taken together, these validations demonstrated a complete N. benthamiana reference genome with the highest quality of all versions.

Genomic landscape of the complete centromeres

Centromeres usually contain megabase-scale repetitive sequences to support CENH3 occupancy and loading⁶. Here we determined the N. benthamiana centromeres based on CENH3 enrichment via aligning the CENH3 ChIP-seq reads to NB.T2T assembly (Supplementary Fig. 7 and Supplementary Table 6). Most identified centromeres showed one dominant peak of CENH3 enrichment, whereas CEN01, CEN10 and CEN14 showed two or three adjacent peaks (Extended Data Fig. 4). The centromere length varied from 0.30 Mb (CEN06) to 5.04 Mb (CEN04) with an average size of 2.03 Mb, and the average GC content varied between 35.11% (CEN02) and 57.73% (CEN16). Based on the ratio of long and short arms, we determined 14 acrocentric (t), two subtelocentric (st) and three submetacentric (sm) chromosomes in N. benthamiana (Fig. 3a,b).

To validate centromere assembly, we aligned ONT and HiFi reads to the NB.T2T assembly, which showed uniform coverage across the centromeres (Fig. 3c). We also examined our Bionano optical mapping across the assembly, of which the optical contigs were consistent with the structure of satellite-free centromeres (Supplementary Fig. 8), while the satellite-enriched centromeres were not fully covered by the optical contigs due to lack of labeling sites in the satellite arrays (Supplementary Fig. 9). The 19 centromeres were mainly classified into three distinct types: satellite-type, Gypsy-type, and NUMTs-type, based on their composition of repetitive elements (Supplementary Table 7). The four satellite-type centromeres (CEN04, CEN08, CEN16 and CEN17) contained long arrays of satellite repeat, resulting in higher GC content and longer centromere size, whereas the seven Gypsy-type centromeres (CEN06, CEN07, CEN11, CEN12, CEN13, CEN15 and CEN19) were dominated by Ty3/Gypsy retrotransposons with lower GC content and shorter length (Fig. 3b). Noteworthy, the combined-type centromeres (CEN01, CEN05, CEN10 and CEN14) always presented two primary CENH3 peaks that were enriched with satellites and retrotransposons, respectively, while another three chromosomes showed merged pattern of both satellite- and Gypsy-type centromeres (CEN03, CEN09 and CEN18), therefore exhibiting comparable levels of satellites and Ty3/Gypsy retrotransposons. The last centromere of CEN02 was mainly composed of LTR/Caulimovirus (37.85%) and NUMTs (27.10%), which was identified as NUMTs-type as the CENH3 mainly enriched on the NUMT sequences (Extended Data Fig. 5a). The presence of mixed centromeres within a single genome probably indicates the centromeric sequences are under rapid and ongoing evolution in N. benthamiana, which will facilitate the study on evolutionary context of these mechanisms in the future.

Origin of the centromeric satellites in N. benthamiana

Despite highly conserved centromeric function in all eukaryotes, centromere sequences show remarkable diversity across species³. To understand this centromere evolution, we dissected the centromeric composition from multiple representative species, nine from T2T assemblies and four from high-quality draft assemblies with a phylogenetic tree built using their CENH3 proteins (Fig. 4a and Supplementary Fig. 10). The centromeres in Arabidopsis⁶, rice⁵ and soybean²⁴ have long satellite arrays. By contrast, Solanaceae centromeres are dominated by Ty3/Gypsy retrotransposons^11–14, similar to a closely-related species of Coffea arabica²⁵, suggesting the Gypsy-rich centromere might be a common feature of Solanaceae species (Fig. 4a). Surprisingly, unlike other Solanaceae members, N. benthamiana not only included Gypsy-type centromeres, but also contained multiple satellite-type centromeres, with satellite repeats accounting for 55.3% of whole centromere DNA (Fig. 4a). We also identified low copies of the same centromeric satellites in the N. benthamiana wild accession (QLD) (Supplementary Fig. 11), indicating these centromere-specific satellites emerged prior to the divergence between LAB and QLD accessions at around 1.1 Mya²⁶.

We further annotated genome-wide tandem repeats using TRASH³, and found 11 centromeres were mainly featured by 33-bp and 43-bp microsatellites, in addition to a small amount of 209-bp and 448-bp satellite (Fig. 4b; Extended Data Fig. 6a; Supplementary Fig. 12). These satellites evolved through tandem duplication leading to higher-order repeats as 76-bp, 86-bp, 162-bp and 657-bp monomer. The centromeric satellite arrays spanned on average 2.20 Mb, ranging from 0.56 Mb (CEN18) to 5.0 Mb (CEN04) per chromosome. Interestingly, we found these centromeric satellites were related to the 45S rDNA and 5S rDNA sequences (Fig. 4c). For instance, the 33-bp and 43-bp monomer share a common motif with the subrepeats of 25-18S rDNA intergenic spacer (IGS), emerging in all centromeric satellite arrays with a total length of 21.67 Mb. The 209-bp and 448-bp satellite are identical to whole unit of 5S rDNA with different IGS sequences, mainly located in centromeric and centromere-proximal regions (Supplementary Fig. 13). Yang et al. also reported that a 5S rDNA array was recruited as functional centromere in one chromosome of the tetraploid switchgrass²⁷. Collectively, these results indicated that the 45S rDNA subrepeats and the 5S rDNA arrays could be recruited as centromeric satellites in N. benthamiana.

Centromere invaded by distinct retrotransposons

Although Ty3/Gypsy accounted for 42.4% of the genome, only 597 elements were identified as complete, belonging to seven subfamilies as Athila, CRM, Galadriel, Ogre, Reina, Retand and Tekay (Extended Data Fig. 6b). Among which, 75 retrotransposons were located within the centromeres, particularly within the satellite arrays, including 8 CRM and 67 Galadriel with a mean length of 7.04 and 5.79 kb, respectively (Supplementary Table 8). Phylogenetics of 279 genome-wide Galadriel retrotransposons showed five distinctive groups, with centromeric elements mainly belonging to younger groups of Galadriel2B, Galadriel3 and Galadriel4 (Supplementary Fig. 14). Analyzing the insertion time of centromeric Galadriel elements showed that at least three rounds (~ 0.5, ~ 1.5, ~ 3.0 Mya) of extensive invasions were occurred during the burst of LTR/Gypsy retrotransposons in whole genome (Supplementary Fig. 15). The recent round of invasion accounted for 65% of intact Galadriel in centromeres, a common feature for both LAB and QLD accessions. The centromere-specific retrotransposons have also been identified in other plants. For example, the CRM retrotransposons were enriched in maize¹⁰, potato¹¹ and pepper¹⁴ centromeres, whereas Athila and Tekay retroelements coexisted with CENH3 in Arabidopsis⁶ and tobacco¹³ centromeres, respectively. These facts imply that the invaded LTR retrotransposons contribute to centromere function, although the mechanisms remain unknown.

Considering that the intact LTR-RTs only contribute to a small portion of total retrotransposons, we also investigated the truncated Gypsy elements in centromeric regions (Supplementary Table 9 and Supplementary Fig. 16). The satellite-type centromeres were mainly invaded by CRM and Galadriel, while Gypsy-type centromeres were predominately occupied by Ogre and Tekay. The combined- and merged-type centromeres showed a mixture of these four subfamilies. These findings imply that the CRM/Galadriel and Ogre/Tekay retrotransposons specifically invaded satellite-type and Gypsy-type centromeres, respectively, shaping their centromere architectures and driving their divergent evolution in tetraploid N. benthamiana.

NUMTs surround the maternal centromeres

Chloroplast (cpDNA) and mitochondrial DNA (mtDNA) are frequently transferred into the nuclear genome during the genome evolution, which still seems to be an ongoing process²⁸. The nuclear transfer of organelle DNA, accumulated mainly at heterochromatin such as the pericentromeric and subtelomeric regions, acting as a main evolutionary force driving genomic instability and diversity²⁸. We first assembled the chloroplast and mitochondrial genomes (Supplementary Fig. 17), searched the foreign organelle fragments in NB.T2T nuclear genome, and found more abundant mtDNA than cpDNA at peri/centromeric regions (Supplementary Fig. 13 and Supplementary Table 10). Totally, 113 inserted NUMTs in length of 7.44 Mb were identified genome-wide, among which 81 segments were located at (0.85 Mb) or close to (5.03 Mb) the centromeres, suggesting the 5-Mb regions proximity to the centromeres are hotspots for mtDNA integration. The frequent colocalization of NUMTs at peri/centromeres was probably due to the stable heterochromatic environment facilitating the integration and elimination of the organelle fragments^28,29.

The 19 chromosomes of allotetraploid N. benthamiana were later phased into A and B subgenomes (see below). Interestingly, 9 out of 10 chromosomes in B subgenome had mtDNA integration at peri/centromeric regions, which only occurred in 2 out of 9 centromeres in A subgenome, demonstrating an obvious subgenome bias of mtDNA integration in N. benthamiana (Supplementary Tables 7 and 10). Given the maternal inheritance of mitochondria, we propose the B subgenome might be maternal derived, agreeing with a previous report³⁰. In addition, the centromeric NUMTs were frequently identified at the boundary of long satellite arrays (Supplementary Fig. 13). For centromere CEN02, CENH3 was exclusively enriched on the NUMTs, but dropped sharply at arrays of LTR/Caulimovirus (Extended Data Fig. 5a), suggesting the foreign mtDNA insertions acting as functional centromere in N. benthamiana. Interestingly, the six NUMTs insertions on CEN08 were evenly distributed (Extended Data Fig. 5b), probably due to the HORs and NUMTs formed “composite satellites”, which then expanded through tandem duplications, resembling the multiple large duplications of AtCEN178 arrays in A. thaliana³. The four NUMTs of the quadruple duplication had strong collinearity due to a common origin, although surprisingly CENH3 ChIP-seq signals peaked at the first two NUMTs but dropped at the latter two for unknown reason. Together, these results revealed the NUMTs played vital roles in diversifying and shaping the centromere architecture in N. benthamiana.

Precise kinetochore assembly sites in N. benthamiana

The functional kinetochore assembly on centromere is essential to attach chromosomes to spindle microtubules in eukaryotes, facilitating the chromosome segregation during mitosis and meiosis¹. As CENH3 defines the functional centromere, we investigated how CENH3 enrichment related to elements as centromeric satellites, Ty3/Gypsy retroelements and NUMTs (Extended Data Fig. 6c). By sorting the 10-kb bins of centromere sequences based on their CENH3 enrichment levels measured as log2(ChIP/control), we observed that majority of satellite arrays showed a medium CENH3 enrichment (0 ~ 2). High-enrichment regions (> 3) were dominated by Gypsy retrotransposons (44.8%), and then followed by satellite repeats (25.8%) and minor NUMTs (13.6%), suggesting the differential involvement of all three elements in centromere function. By contrast, CENH3 enrichment was significantly lower within LTR/Caulimovirus (< 0) than other elements, indicating that the LTR/Caulimovirus invasion is detrimental to centromere function in N. benthamiana. Overall, CENH3 signals were more likely colocalized with Gypsy retrotransposons rather than satellite arrays, although satellite repeats occupied 55.3% of centromeres (Fig. 4a). Fine-scale profiling of ChIP-seq around feature sets confirmed the highest enrichment for CRM and Tekay retrotransposons, indicating that the CRM and Tekay elements were the precise CENH3 binding sites and potentially shaped the centromeric architecture in satellite-type and Gypsy-type centromeres, respectively (Extended Data Fig. 6d). Our finding was contrary to the report in Arabidopsis⁶, where CENH3 enrichment was exclusively associated with CEN180 arrays, and drastically dropped on Athila retrotransposons, suggesting distinctive roles of centromeric satellites and retrotransposons in Arabidopsis and N. benthamiana.

Centromeres contain a characteristic combination of epigenetic marks that are critical for its function. To assess the epigenetic features of the centromeres in N. benthamiana, we re-analyzed the recently published histone ChIP-seq dataset¹⁵. As expected, the transcriptionally active mark of H3K4me3 was closely correlated with gene density, while the heterochromatin marker of H3K9me2 was enriched on genomic repeat, as previously reported¹⁵. However, the centromeres of N. benthamiana displayed distinct patterns of histone methylation marks compared to Arabidopsis⁶. The Arabidopsis centromeric arrays were enriched in H3K9me2, and depleted in H3K4me3 and H3K27me3, whereas the centromeric satellites here in N. benthamiana were enriched in H3K4me3 and H3K9me2, and depleted in H3K27me3 (Fig. 3c and Extended Data Fig. 6b). Furthermore, the euchromatin specific active marker of H3K4me3 here was enriched on centromeric satellites, while the heterochromatin specific repressive H3K9me2 here was depleted on centromeric transposons, although the underlying mechanism is unknown (Extended Data Fig. 6d). Other centromeric components as NUMTs and LTR retrotransposons showed no obvious enrichment of histone modifications. The histone ChIP-seq analysis demonstrated that the H3K4me3, H3K9me2 and H3K27me3 were all depleted on centromeric NUMTs, exhibiting an epigenetic pattern differing from that of centromeric satellites (Extended Data Fig. 5a,b). The distinctive histone modification patterns in three types of centromeres indicated that they had unique epigenetic regulations in N. benthamiana.

A model for centromere formation and evolution in N. benthamiana

Centromeres are composed of long satellite arrays in many plant and animal species such as rice⁵, Arabidopsis⁶ and humans⁴, while the Solanaceae plants such as potato¹¹, tomato¹² and pepper¹⁴ are featured by Ty3/Gypsy dominated centromeres. The emergence of satellite centromeres has been elusive, probably evolves from “neocentromeres” that originally do not contain satellite repeats². Hence, the N. benthamiana provides a model system for studying centromere evolution as it possessed satellite-based and satellite-free centromeres.

Within the N. benthamiana tetraploid genome, we searched for syntenic blocks that embedded whole centromere, and found two cenhaps (centromere-spanning haplotype blocks)³ across CEN01 (Fig. 4d,e) and CEN10 (Extended Data Fig. 7a,b). These two centromere-embedded syntenic blocks were also conserved among tobacco (SubT and SubS genome), potato and pepper. However, the centromeres of these species were outside the syntenic blocks, suggesting the centromeres of CEN01 and CEN10 in N. benthamiana were born recently after diverged with tobacco. Thus, based on the findings from centromere analysis, we propose a model for the evolution of satellite-based centromeres in N. benthamiana which involves two processes: the formation of neocentromeres at non-centromeric regions by acquiring CENH3 protein, and then transition to mature centromeres by satellite expansion, centromeric retrotransposon integration, and sometimes mtDNA invasion.

Subgenome assignment of the allotetraploid N. benthamiana

N. benthamiana is an allotetraploid plant indigenous to Australia, belonging to the Section Suaveolentes of Nicotiana, whose subgenome assignment has been contentious. Although most studies propose Noctiflorae/Petunioides and Alatae/Sylvestres sections to be their diploid ancestors, chromosomal assignments for two subgenomes have showed conflicting results^15,30,31. To solve this controversy, we firstly phased subgenomes of the allotetraploid N. benthamiana by using k-mer-based method. We used the SubPhaser³² to identify subgenome-specific signatures (Extended Data Fig. 8a,b and Supplementary Table 11), and found 51,610 differential enriched k-mers, which confidently divided the chromosomes into two distinct groups (Fig. 5a). Meanwhile, principal component analysis (PCA) separated nine and ten chromosomes into A (SubA) and B (SubB) subgenome, respectively, with PC1 mainly explaining 56.3% of the variance between the two subgenomes (Fig. 5b). The subgenomic assignments here conflicted with a recently published study¹⁵, with only 12 of the 19 chromosomes were partitioned consistently (Fig. 5c). As the previously used disjoint subset partitioning method is not reliable, we suggested the results here generated by k-mer-based approach are more representative of the subgenomic assignment in N. benthamiana. But it is worth noting that the chromosomes in Section Suaveolentes are formed by a mosaic of both subgenomes³³, of which some chromosomes are “enriched” in one or another parent (Extended Data Fig. 8a,b), so the subgenome assignment in N. benthamiana is only a simplified strategy.

We then used phylogeny-based and read mapping methods to infer the potential diploid ancestors. Phylogenetic relationships were constructed using a concatenation of 2,840 single-copy orthologues with potato and tomato as outgroup (Fig. 5d). Our data supported the clustering of A subgenome with N. sylvestris, whereas B subgenome was sister to N. glauca and N. attenuata. We then assessed the genomic contributions of diploid species to the allotetraploid by mapping Illumina data of N. benthamiana to six Nicotiana diploid genomes (Fig. 5e). The read depth in three diploid species of N. sylvestris, N. glauca and N. attenuata was higher than that in N. longiflora and N. tomentosiformis. The mapping of Illumina reads of the N. sylvestris and N. attenuata to the allotetraploid N. benthamiana genome also confirmed the A and B subgenome as N. sylvestris and N. attenuata related (Fig. 5f). Synteny analysis between the two subgenomes of N. benthamiana and two ancestry relatives of N. sylvestris and N. attenuata revealed remarkable synteny blocks, although extensive recombination events had occurred between and within two subgenomes (Fig. 5g and Extended Data Fig. 8c,d). The BUSCO score of two subgenomes only reached ~ 80% of the complete genes (Supplementary Fig. 18), suggesting the allotetraploid lost many homeologous sequences during the process of “genomic shock”.

Paleohistory of the N. benthamiana genome

Based on our subgenome analyses and previous reports^{20,30,31,34,35}, we summarized the paleohistory of the N. benthamiana genome evolution (Fig. 6 and Supplementary Fig. 19). As members of the nightshade family, the Nicotiana and Solanum species shared a common whole genome triplication (WGT) event. Then the two genera diverged approximately at 23–44 Mya (http://www.timetree.org/), consistent with our phylogenomic analysis of Nicotiana spp. and related angiosperm species (Supplementary Fig. 19). We phased the two subgenomes of N. benthamiana, and inferred the potential paternal and maternal ancestors as sisters to N. sylvestris and N. attenuata/N. glauca, respectively, corresponding to Sylvestres and Noctiflorae/Petunioides clades. Due to lacking chromosome-level genome of N. glauca, only N. attenuata was used for subsequent syntenic analysis and demonstration. The divergence between Sylvestres and Noctiflorae/Petunioides was estimated to be 6–10 Mya^31,34,35. Then the Section Suaveolentes arose from a hybridization event during the Pliocene transition at 5–6 Mya^30,35. Afterwards, the allopolyploid experienced extensive rearrangements between subgenomes, termed as ‘genomic shock’ commonly found in polyploid evolution^35,36. The long-term genome diploidization in allopolyploid Nicotiana eventually triggered chromosome number reduction, genome downsizing, centromere gain and loss, satellite emergence, and etc.³⁷, giving rise to the recent N. benthamiana species (Fig. 6).

The application of Nicotiana benthamiana as a model plant has increased exponentially over the past decade. Leveraging advances in PacBio HiFi and ultra-long ONT sequencing technologies, we have achieved a complete T2T assembly of all 19 chromosomes of N. benthamiana, providing a valuable resource for plant research community. We characterized three distinct types of centromeres in N. benthamiana through CENH3 ChIP-seq profiling, and resolved the epigenetic landscape within the centromeres using public histone ChIP-seq data. Contrary to well-studied Arabidopsis and human centromeres, the N. benthamiana centromeres showed different genetic and epigenetic patterns. In humans, the centromeres are dominated by 171-bp α-satellite repeat, which can be organized in long arrays varying from 340 kb to 4.8 Mb⁴. And the strongest CENP-A binding is within active α-satellite arrays with low CpG methylation⁴. The Arabidopsis centromeres contain megabase-scale 178-bp satellite repeats, characterized by densely DNA methylation, H3K9me2 enrichment and H3K4me3 depletion⁶. By contrast, our study of N. benthamiana revealed its satellite-type centromeres were featured by 33-/43-bp microsatellites with enriched H3K9me2 and H3K4me3 marks, while the Gypsy-type centromeres were dominated by Ty3/Gypsy retrotransposons without obvious epigenetic marks. Both types of centromeres showed high CpG methylation levels. These results demonstrate centromere diversity on genetic and epigenetic landscapes among species.

Satellite repeats have been a characteristic of centromeres in many plant and animal species. However, the satellite-free centromeres are commonly found in Solanaceae plants such as potato¹¹, tomato¹², tobacco¹³ and pepper¹⁴. We previously found pepper centromeres were enriched with Gypsy especially CRM retrotransposons, and their abundance varied among different Capsicum species¹⁴. The recent potato gap-free genome also showed that its centromeres contained 49.2% of Gypsy LTR-RTs¹¹, among which six centromeres contained long arrays of satellite repeats that mostly originated from Gypsy retrotransposons³⁸. The repeat of Sv98 in potato produced strong FISH (Fluorescence in situ hybridization) signals in most centromeres, which was further defined as a centromere-targeted CRM retrotransposon³⁸. Similarly, the Gypsy retrotransposons of TGR4 and Nt2-7/Nto1 were exclusively found in the tomato¹² and tobacco¹³ centromeres, respectively, which were annotated as Tekay elements. It seems that the CRM or Tekay elements probably shape the architecture of Solanaceae centromeres, acting as functional components of centromeres. This is consistent with our finding that the centromeric CRM and Tekay retrotransposons are preferred sites of CENH3 binding in N. benthamiana, although the mechanisms of the retrotransposon functioning as centromere remains unknown.

Centromeric satellites with unit ranging from 150 bp to 180 bp are highly phased with the CENH3 nucleosomes, as a monomer of satellite is desirable for wrapping a single nucleosome^1,27. Interestingly, the lengths of satellite monomer in N. benthamiana are out of this range, being as low as 33-/43-bp or as high as 209-/448-bp, which were only restricted to 11 centromeres and have not spread to the other eight centromeres. Although a few studies reported satellite repeats in Solanaceae plants of potato³⁸ and tomato¹², we regarded it here as enrichment of centromeric retrotransposons rather than the canonical satellite repeats. We proposed the centromeric 33-/43-bp satellite repeats newly emerged in N. benthamiana, as these arrays were absent in other closely related species. As three distinct types of centromeres are present in N. benthamiana, an intriguing question is whether the centromere evolution involves transition from one state to another. Jiang and colleagues hypothesized that young neocentromeres were satellite-free, then invaded by satellite repeats during evolution, and eventually resulted in mature satellite-based centromeres^1,2. If this was the case, the neocentromeres of N. benthamiana were most likely formed firstly during tetraploidization process, as frequent fissions and fusions had occurred within two ancestral subgenomes (Fig. 6). The fused chromosomes need a functional centromere to assemble kinetochore and attach spindle fiber¹. Afterwards, the neocentromeres diverged with the expansion of distinct retrotransposons. Then satellite repeats emerged and amplified exclusively in CRM/Galadriel enriched centromeres. We speculated the invasions of centromeric satellites into the centromeric retrotransposons rather than the opposite situation, considering that the centromeric retrotransposons served as the core of centromeres and the satellite repeats evolved more rapidly than retrotransposons. Why the satellite repeats did not invade Tekay/Ogre enriched centromeres needs further investigation in the future.

The presence of mixed or multiple types of centromeres within a single genome has also been reported in rice³⁹ and maize¹⁰. The merged-type centromeres here are proposed to represent an intermediate stage of the centromeric satellite expansion between the satellite-free neocentromeres and mature satellite-based centromeres. A mechanism behind this is the ‘female meiotic drive’, of which the emergence of new satellite array or expansion of an existing centromeric satellite would recruit more kinetochore proteins and confer an advantage during chromosome segregation in female meiosis^40,41. The DNA transfer from organelles to nucleus has been reported in many species, which were frequently inserted at the pericentromeric regions of human⁴², Arabidopsis⁴³ and rice²⁹. Puertas and González-Sánchez explained this phenomenon as the repair of double-strand breaks via the non-homologous end-joining mechanism that incorporated the mtDNA into the nuclear genome²⁸. Previous studies on human⁴ and Arabidopsis⁶ centromeres lacked the characterizations of the genomic and epigenomic patterns of NUMTs in the vicinity of centromeres. The study here not only identified frequent NUMT insertions surrounding the maternal-derived centromeres, but also found centromeric NUMTs evolving to the core of functional centromere. The novelty of this finding lies in enriching our knowledge on effects of foreign mtDNA on nuclear genome. Taken together, the genetic and epigenetic landscape of centromeres will facilitate the understanding of centromere and polyploid genome evolution in Solanaceae and other eukaryotes.

Acknowledgments

We would like to thank the Bioinformatics Platform at Peking University Institute of Advanced Agricultural Sciences for providing the high-performance computing resources. This work was supported by the Key R&D Program of Shandong Province (ZR202211070163), Taishan Scholars Program and Natural Science Foundation for Distinguished Young Scholars (ZR2023JQ010) of Shandong Province.

Author contributions

L.G. conceived and supervised the project. W.C. performed genome assembly and analysis of centromeric sequences. M.Y. prepared figures and tables. S.C., J.S. and J.W. performed epigenetic analysis. D.M. conducted epigenome sequencing. J.L. and Z.Z. generated genome sequencing data. W.C. and L.G. wrote the manuscript. All authors read and approved the final version of the manuscript.

Competing interests

The authors declare no competing interests.

Data availability

The raw sequencing data and the genome assembly have been deposited at the China National Center for Bioinformation (https://ngdc.cncb.ac.cn/) under project number PRJCA022857. The nuclear genome assembly and annotation are also available at https://github.com/Weikai-47/Nicotiana_Project. The mitochondria and chloroplast genomes have been deposited at the China National Center for Bioinformation under accession number C_AA066595 and C_AA066594, respectively.

Code availability

This manuscript does not report original code.

Jiang, J., Birchler, J.A., Parrott, W. A. & Dawe, R. K. A molecular view of plant centromeres. Trends Plant Sci.8, 570-575 (2003).
Zhang, H. et al. Boom-bust turnovers of megabase-sized centromeric DNA in Solanum species: rapid evolution of DNA sequences associated with centromeres. Plant Cell26, 1436-1447 (2014).
Wlodzimierz, P. et al. Cycles of satellite and transposon evolution in Arabidopsis centromeres. Nature618, 557-565 (2023).
Altemose, N. et al. Complete genomic and epigenetic maps of human centromeres. Science376, eabl4178 (2022).
Shang, L. et al. A complete assembly of the rice Nipponbare reference genome. Mol. Plant16, 1232-1236 (2023).
Naish, M. et al. The genetic and epigenetic landscape of the Arabidopsis centromeres. Science374, eabi7489 (2021).
Zhao, J. et al. Centromere repositioning and shifts in wheat evolution. Plant Commun. 4, 100556 (2023).
Ahmed, H. I. et al. Einkorn genomics sheds light on history of the oldest domesticated wheat. Nature620, 830-838 (2023).
Liu, Q. et al. Non–B-form DNA tends to form in centromeric regions and has undergone changes in polyploid oat subgenomes. Proc. Natl. Acad. Sci. USA120, e2211683120 (2023).
Chen, J. et al. A complete telomere-to-telomere assembly of the maize genome. Nat. Genet.55, 1221-1231 (2023).
Yang, X. et al. The gap-free potato genome assembly reveals large tandem gene clusters of agronomical importance in highly repeated genomic regions. Mol. Plant16, 314-317 (2023).
Chang, S. B. et al. FISH mapping and molecular organization of the major repetitive sequences of tomato. Chromosome Res.16, 919-933 (2008).
Nagaki, K. et al. Coexistence of NtCENH3 and two retrotransposons in tobacco centromeres. Chromosome Res.19, 591-605 (2011).
Chen, W. et al. Two telomere-to-telomere gapless genomes reveal insights into Capsicum evolution and capsaicinoid biosynthesis. Nat. Commun.15, 4295 (2024).
Ranawaka, B. et al. A multi-omic Nicotiana benthamiana resource for fundamental research and biotechnology. Nat. Plants9, 1558-1571 (2023).
Cheng, H., Concepcion, G. T., Feng, X., Zhang, H. & Li, H. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nat. Methods18, 170-175 (2021).
Bombarely, A. et al. A draft genome sequence of Nicotiana benthamiana to enhance molecular plant-microbe biology research. Mol. Plant Microbe Interact.25, 1523-1530 (2012).
Kurotani, K. I. et al. Genome sequence and analysis of Nicotiana benthamiana, the model plant for interactions between organisms. Plant Cell Physiol.64, 248-257 (2023).
Wu, Y. et al. Phylogenomic discovery of deleterious mutations facilitates hybrid potato breeding. Cell186, 2313-2328 (2023).
Wang, J. et al. High-quality assembled and annotated genomes of Nicotiana tabacum and Nicotiana benthamiana reveal chromosome evolution and changes in defense arsenals. Mol. Plant17, 423-437 (2024).
Rabanal, F. A. et al. Pushing the limits of HiFi assemblies reveals centromere diversity between two Arabidopsis thaliana genomes. Nucleic Acids Res.50, 12309-12327 (2022).
Rhie, A., Walenz, B. P., Koren, S. & Phillippy, A. M. Merqury: reference-free quality, completeness, and phasing assessment for genome assemblies. Genome Biol.21, 245 (2020).
Ou, S., Chen, J. & Jiang, N. Assessing genome assembly quality using the LTR Assembly Index (LAI). Nucleic Acids Res.46, e126 (2018).
Wang, L. et al. A telomere-to-telomere gap-free assembly of soybean genome. Mol. Plant16, 1711-1714 (2023).
de Castro Nunes, R. et al. Structure and distribution of centromeric retrotransposons at diploid and allotetraploid Coffea centromeric and pericentromeric regions. Front. Plant Sci.9, 175 (2018).
Cauz-Santos, L. A. et al. Genomic insights into recent species divergence in Nicotiana benthamiana and natural variation in Rdr1 gene controlling viral susceptibility. Plant J.111, 7-18 (2022).
Yang, X. et al. Amplification and adaptation of centromeric repeats in polyploid switchgrass species. New Phytol.218, 1645-1657 (2018).
Puertas, M. J., & González-Sánchez, M. Insertions of mitochondrial DNA into the nucleus-effects and role in cell evolution. Genome63, 365-374 (2020).
Matsuo, M., Ito, Y., Yamauchi, R. & Obokata, J. The rice nuclear genome continuously integrates, shuffles, and eliminates the chloroplast genome to cause chloroplast-nuclear DNA flux. Plant Cell17, 665-675 (2005).
Schiavinato, M., Marcet‐Houben, M., Dohm, J. C., Gabaldón, T., & Himmelbauer, H. Parental origin of the allotetraploid tobacco Nicotiana benthamiana. Plant J.102, 541-554 (2020).
D’Andrea, L. et al. Polyploid Nicotiana section Suaveolentes originated by hybridization of two ancestral Nicotiana clades. Front. Plant Sci. 14, 999887 (2023).
Jia, K. H. et al. SubPhaser: a robust allopolyploid subgenome phasing method based on subgenome-specific k-mers. New Phytol.235, 801-809 (2022).
Chase, M. W. et al. Molecular systematics, GISH and the origin of hybrid taxa in Nicotiana (Solanaceae). Ann. Bot.92, 107-127 (2003).
Wang, S. et al. Phylotranscriptomics supports numerous polyploidization events and phylogenetic relationships in Nicotiana. Front. Plant Sci.14, 1205683 (2023).
Clarkson, J. J., Dodsworth, S. & Chase, M. W. Time-calibrated phylogenetic trees establish a lag between polyploidisation and diversification in Nicotiana (Solanaceae). Plant Syst. Evol.303, 1001-1012 (2017).
Lim, K. Y. et al. Sequence of events leading to near-complete genome turnover in allopolyploid Nicotiana within five million years. New Phytol.175, 756-763 (2007).
Koukalova, B. et al. Fall and rise of satellite repeats in allopolyploids of Nicotiana over c. 5 million years. New Phytol.186, 148-160 (2010).
Gong, Z. et al. Repeatless and repeat-based centromeres in potato: implications for centromere evolution. Plant Cell24, 3559-3574 (2012).
Song, J. et al. Two gap-free reference genomes and a global view of the centromere architecture in rice. Mol. Plant14, 1757-1767 (2021).
Malik, H. S. & Henikoff, S. Major evolutionary transitions in centromere complexity. Cell138, 1067-1082 (2009).
Naish, M. & Henderson, I. R. The structure, function, and evolution of plant centromeres. Genome Res. 34, 161-178 (2024).
Wei, W. et al. Nuclear-embedded mitochondrial DNA sequences in 66,083 human genomes. Nature611, 105-114 (2022).
Michalovová, M., Vyskot, B., & Kejnovsky, E. Analysis of plastid and mitochondrial DNA insertions in the nucleus (NUPTs and NUMTs) of six plant species: size, relative age and chromosomal localization. Heredity111, 314-320 (2013).

Plant materials and genome sequencing

The N. benthamiana LAB strainwas routinely grown in a custom soil mix (potting mix, clay and vermiculite) with a day/night cycle of 16 h/8 h at a constant temperature of 25°C in greenhouse of Peking University Institute of Advanced Agricultural Sciences, Weifang, China. After two weeks, the fresh young leaves were collected and used for genome sequencing and optical mapping. Tissues of roots, leaves, stems, flowers and seeds at two days post anthesis were used to perform Illumina RNA-sequencing on Novaseq 6000 platform and full-length transcriptome sequencing on PacBio Sequal IIe instrument. The high-molecular-weight (HMW) genomic DNA were extracted and used for the construction of sequencing library⁴⁴. To achieve telomere-to-telomere assembly of N. benthamiana, we performed 110× PacBio HiFi sequencing on Sequel II platform, 40× ONT ultralong sequencing on GridION X5/PromethION sequencer, and 100× Illumina paired-end sequencing on Novoseq 6000 platform. Then we constructed Hi-C sequencing libraries using HindIII restriction enzyme, which were subsequently sequenced on Illumina NovaSeq 6000 to obtain 150× paired-end short reads. The Hi-C data were classified as valid or invalid interactions using HiC-Pro (v3.1.049)⁴⁵, and only valid interaction pairs were retained for subsequent analysis. Details of genome and transcriptome sequencing are described in the Supplementary Methods.

Bionano optical map generation

The fresh young leaves were collected from N. benthamiana and processed in Grand Omics (Wuhan, China) for construction of optical maps. The high-molecular-weight DNA was extracted using the Prep Plant DNA Isolation Kit (Bionano Genomics), and labelled using the DLE-1 enzyme following the DLS protocol (Bionano Genomics). The labeled DNA was then loaded onto the Saphyr system (Bionano Genomics) for the imaging analysis. A total of 329.6× coverage of effective molecules was generated with an average label density of 14.96/100 kb. The DNA molecules were subsequently assembled using Bionano Solve (v3.7) and the Bionano images of the assembly were visualized using Bionano Access (v1.7.1). Details of Bionano data processing are described in the Supplementary Methods.

Genome assembly

To obtain a high contiguity assembly, Hifiasm (v0.19.5)¹⁶ was employed combining the HiFi and ultralong ONT reads with the parameters of “-l 0 --n-hap 2 --hom-cov=111”. The well-assembled HiFi & ONT assembly was in a total length of 2.85 Gb with contig N50 of 144.7 Mb, leaving only three gaps. Two gaps located in 45S rDNA repeat arrays of Chr06 and Chr19 were resolved by following NOR assembly (Supplementary Methods). The last gap, corresponding to NUPT insertion in Chr12, was closed by ultralong ONT reads and confirmed by coverage depth in the IGv⁴⁶. After gap filling, we generated an assembly comprising 19 contigs with 38 telomeres. Then Hi-C reads were used to anchor all contigs using the pipeline of BWA (v0.7.17)⁴⁷, Juicer (v1.5)⁴⁸ and 3D-DNA (v180419)⁴⁹. We manually checked the assembly through chromatin interaction patterns in Juicebox (v1.11.08)⁵⁰. An mis-assembly in Chr03 was corrected and the adjustment was then confirmed in the IGv⁴⁶. Finally, we obtained a T2T genome assembly of N. benthamiana with total length of 2.85 Gb and contig N50 of 146.4 Mb.

Validation of the T2T genome assembly

To assess the quality of the genome assembly, we performed extensive validation in multiple approaches. Firstly, Bionano optical map were aligned to the in silico genome assembly using RefAlinger (v12432.12542rel) tool in the pipeline of Bionano Solve (v3.7). Then the Hi-C interaction matrix of the final assembly was manually checked in Juicebox⁵⁰ (Supplementary Methods). Secondly, we aligned the NB.T2T assembly with recently published assemblies of N. benthamiana using Syri (v1.6.3)⁵¹ and found good collinearities. Thirdly, we mapped HiFi, ONT and NGS reads against the NB.T2T assembly and checked for read coverage for any abnormal signals. Finally, to assess genome completeness, we applied BUSCO (v5.4.3)⁵² for ortholog detection using solanales_odb10 database (n=5,950). Quality value (QV) was estimated using Merqury (v1.3)²² from HiFi reads. LTR Assembly Index (LAI) was calculated using LTR_retriever (v2.9.0)⁵³ to evaluate the assembly continuity of LTR retrotransposons²³.

Repeat sequence annotation

The de novo repeat library for the N. benthamiana was constructed by RepeatModeler (https://github.com/Dfam-consortium/RepeatModeler). Then combining the universal Repbase database (version 20181026), the repetitive elements were annotated using RepeatMasker (v4.1.2)⁵⁴. The intact LTR elements were identified using the pipeline of LTR_Finder (v1.2)⁵⁵, LTRharvest (v1.6.2)⁵⁶, and LTR_retriever (v2.9.0)⁵³. Then LTR sequences from RepeatMasker and intact LTR elements from LTR_retriever were transferred to TEsorter (v1.3)⁵⁷ to classify the subfamily of LTR-RTs. Among which, the Ty3-Gypsy elements was classified into seven clades as Athila, CRM, Galadriel, Ogre, Reina, Retand, and Tekay. The insertion time of intact LTR retrotransposons were calculated using LTR_retriever⁵³. The TRASH³ pipeline was used to identify satellite repeats. Two major microsatellites were identified in centromeres, including CEN33 (ACGAGTCAGG ACGTGGCAGG ACATGGCCATGGC) and CEN43 (ACGTGTCAGG ACGCGTCAGG ACGCGTCAGG ACATGGCCATGGC).

Gene model annotation

To predict coding-gene models, we applied MAKER (v2.31.11)⁵⁸ pipeline combining evidence from homology protein, transcript evidence, and ab initio prediction. The proteins used for homology-based prediction were from N. attenuata⁵⁹, N. benthamiana^18,19, S. lycopersicum, S. tuberosum¹¹, and universal Swiss-Prot proteins. The redundant sequences were excluded using cd-hit (v 4.8.1)⁶⁰. RNA-seq reads were mapped to the assembled genome using HISAT2 (v2.2.1)⁶¹, followed by genome-guided transcript assembly by StringTie (v1.13)⁶². The full-length transcriptome reads were processed using the SMRT Analysis software Isoseq3 (https://github.com/PacificBiosciences/IsoSeq). The SNAP⁶³, GeneMark-ET⁶⁴ and AUGUSTUS⁶⁵ models were trained using MAKER2⁵⁸ and BRAKER2 (v2.1.6)⁶⁶ pipeline as previously reported. Finally, the trained models, protein and transcript evidences were integrated into MAKER2 to predict credible gene structures. Gene models were then manually corrected in IGV-GSAman (https://gitee.com/CJchen/IGV-sRNA) with the support of transcript coverage and previous annotations^18,19.

Gene annotation assessment

To provide a valuable resource for plant research community, we compared our genome annotation (NB.T2T) with two recently published versions, NB.PCP¹⁸ and NB.MP²⁰. First, OMArk (v2.0.3)⁶⁷ was performed to assess not only the completeness but also the consistency of the gene repertoire in three versions. The recommended LUCA.h5 database was used to exploit orthology relationships. Second, OrthoFinder (v2.5.4)⁶⁸ was conducted to identify orthogroups using non-redundant protein sequences from three versions.

Organelle genomes assembly

The mitochondria and chloroplast genomes were assembled with ONT data using NextDenovo (v2.5.0)⁶⁹ incorporated in the PMAT pipeline⁷⁰. Bandage software (v0.8.1)⁷¹ was employed to visualize the organelle genomes, and remove any linear organelle fragments and nuclear sequences. Then ONT ultra-long reads containing mitochondrial genes were extracted and subsequently aligned to the mitochondrial genome using blastn within Bandage to resolve repeat regions within the graphical mitochondrial genome. Finally, we obtained the circular mitochondria and chloroplast genomes, annotated using GeSeq (v2.0.3)⁷² online website, and deposited at the China National Center for Bioinformation under accession number C_AA066595 and C_AA066594, respectively. To identify organelle genome insertions within nucleus genomes, blastn searches were performed and the positions were visualized by R package RIdeogram⁷³.

CENH3 ChIP-seq

Previously, based on the amino acid sequence of NbCENH3, two peptides corresponding to the N-terminus (H2N-ARTKHLALRKQSRPPSRPTA-COOH) or full sequence of NbCENH3 was synthesized to produce anti-CENH3 antibodies in rabbit. The ChIP cloning experiment was then conducted on chromatin extract from young leaves using anti-CENH3 antibodies according to a previously reported protocol⁷⁴. The ChIP library was finally amplified with the VAHTS® Universal DNA Library Prep Kit for Illumina V3 (Vazyme ND607), and sequenced on the Illumina Novaseq 6000 platform to produce 150-bp paired-end reads. Since the two antibodies exhibited consistent results, so only the antibody against the full sequence of NbCENH3 was used in this study with two biological replicates.

Centromere identification in public genomes

To extrapolate a phylogenetic story about plant centromeres, we downloaded eight published T2T assemblies and four draft assemblies. The centromeres of Arabidopsis⁶, rice⁵, wheat⁷, einkorn⁸, oat⁹, maize¹⁰, potato¹¹, pepper¹⁴, and soybean²⁴ were identified by CENH3 ChIP-seq, while the grape⁷⁵ and carrot⁷⁶ centromeres were predicted based on distribution of satellite repeats. Additionally, the centromeres of coffee²⁵ and tomato¹² were found by searching centromere-specific retrotransposons as CRC (Centromeric retrotransposons in Coffea) and TGR4 (Ty3-Gypsy type retrotransposons) in whole genomes, respectively. As for tobacco, two centromeric retrotransposons colocalized with NtCENH3 were helped to predicate putative centromere regions¹³. The predicted centromeres of coffee, tomato and tobacco were further confirmed by Hi-C interaction matrixes, as the centromeres always showed obvious inter-chromosome interactions¹⁴.

ChIP-seq data processing

Histone modification ChIP-seq data was downloaded from NCBI database under BioProject PRJNA881799. The CENH3 ChIP-seq were conducted as above mentioned. Paired-end reads were preprocessed with fastp (v0.23.2)⁷⁷ to remove low-quality bases, and then aligned to the NB.T2T assembly using bowtie2 (v2.5.1)⁷⁸ with default settings. Alignments with mapping quality of < 30 were discarded, and all read duplicates were removed using Samtools(v1.10)⁷⁹. For each dataset, the bamCompare tool from Deeptools (v3.5.1)⁸⁰ was used to quantify histone methylation and CENH3 enrichement in bigwig format. The coverage was calculated as the log2 of the ratio (ChIP/control) per 50-bp bin. For profiling of CENH3 occupancy, alignment bedGraphs were also used to calculate log2(ChIP/control) in 10-kb bins. CENH3 domains were identified by comparing the ChIP and input data using MACS2 (v2.2.7.1)⁸¹ with following parameters: “--nomodel --extsize 291 --keep-dup auto -broad”.

ChIP-seq Metaplots

The above bamCompare generated bigwig files were used to calculate the enrichment level across centromeric elements as CEN33/43 satellites, Gypsy retrotransposons, Copia retrotransposons, CENH3 domains, and CRM/Galadriel/Ogre/Tekay elements using computeMatrix tool from Deeptools⁸⁰ in ‘scale-regions’ mode with parameters of “-regionBodyLength 2000 -beforeRegionStartLength 2000 -afterRegionStartLength 2000”. Then metaplots for CENH3, H3K4me3 and H3K9me2 were plotted with plotHeatmap available from Deeptools. The visualization of centromeric and telomeric satellite repeats was accomplished using StainedGlass⁸².

DNA methylation processing

We used Nanopolish (v0.13.2)⁸³ to identify 5-methylcytosine bases in a CpG context in ONT reads. The raw reads were firstly indexed to the fast5 files, and then aligned to the genome using Minimap2⁸⁴ (-ax map-ont). Then Nanopolish was used to call methylation, and the script of calculate_methylation_frequency.py was used to count the methylation frequency at each CpG sites. We further summarized the methylation levels at 10-kb no-overlapping windows.

Subgenome phasing

We used SubPhaser (v1.2.5)³² to phase and partition the subgenomes of allotetraploid N. benthamiana based on differential k-mers among homoeologous chromosome sets. Briefly, the syntenic analysis was firstly performed by JCVI (v1.1.19)⁸⁵ and MCScanX⁸⁶, and genes in syntenic blocks were identified as homeologues. Then SubPhaser searches for the subgenome specific k-mers, and assigns homoeologous chromosomes into two subgenomes. To infer the potential diploid ancestors of N. benthamiana, we applied mapping-based and phylogenetic-based methods. The Illumina short reads of N. benthamiana were mapped to composite reference of five diploid Nicotiana genomes (N. attenuata, N. glauca, N. longiflora, N. sylvestris, N. tomentosiformis) using sppIDer⁸⁷. The mapping rates against each diploid genome indicated the genetic contributions of diploid species to the allotetraploid. Meanwhile, we also mapped the Illumina short reads of N. sylvestris (NCBI BioProject: PRJEB1826) and N. attenuata (NCBI BioProject: PRJNA316810) to NB.T2T assembly. The coverage depth against two subgenomes revealed their evolutionary relationship in a certain degree.

Phylogenetic analysis

OrthoFinder (v2.5.4)⁶⁸ was used to identify orthogroups using non-redundant protein sequences from 21 species (Supplementary Table 12). Each subgenome of N. benthamiana and N. tabacum was considered as an independent pseudospecies. In total, 361 single-copy orthogroups were identified to construct the maximum likelihood tree using RAxML (v8.2.12)⁸⁸. The CodeML and MCMCTree programs in the PAML (v4.9)⁸⁹ were used to analyze amino acid substitution models and estimate divergence times. IQ-TREE (v2.0.3)⁹⁰ was applied to construct the phylogenetic tree of Gypsy retrotransposons based on RT (reverse transcriptase) domains. The maximum likelihood tree of CENH3 proteins was also conducted using IQ-TREE⁹⁰.

Additional references

44. Zhang, M. et al. Preparation of megabase-sized DNA from a variety of organisms using the nuclei method for advanced genomics research. Nat. Protoc. 7, 467-478 (2012).

45. Servant, N. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol. 16, 259 (2015).

46. Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29, 24-26 (2011).

47. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754-1760 (2009).

48. Durand, N. C. et al. Juicer provides a one-click system for analyzing loop-resolution Hi-C experiments. Cell Syst. 3, 95-98 (2016).

49. Dudchenko, O. et al. De novo assembly of the Aedes aegypti genome using Hi-C yields chromosome-length scaffolds. Science 356, 92-95 (2017).

50. Durand, N.C. et al. Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom. Cell Syst. 3, 99-101 (2016).

51. Goel, M., Sun, H., Jiao, W. B. & Schneeberger, K. SyRI: finding genomic rearrangements and local sequence differences from whole-genome assemblies. Genome Biol. 20, 277 (2019).

52. Manni, M., Berkeley, M. R., Seppey, M., Simao, F. A. & Zdobnov, E. M. BUSCO update: novel and streamlined workflows along with broader and deeper phylogenetic coverage for scoring of eukaryotic, prokaryotic, and viral genomes. Mol. Biol. Evol. 38, 4647-4654 (2021).

53. Ou, S. & Jiang, N. LTR_retriever: A highly accurate and sensitive program for identification of long terminal repeat retrotransposons. Plant Physiol. 176, 1410-1422 (2018).

54. Tarailo-Graovac, M. & Chen, N. Using RepeatMasker to identify repetitive elements in genomic sequences. Curr. Protoc. Bioinformatics 5, 4-10 (2009).

55. Xu, Z. & Wang, H. LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 35, W265-W268 (2007).

56. Ellinghaus, D., Kurtz, S. & Willhoeft, U. LTRharvest, an efficient and flexible software for de novo detection of LTR retrotransposons. BMC Bioinformatics 9, 18 (2008).

57. Zhang, R. G. et al. TEsorter: an accurate and fast method to classify LTR-retrotransposons in plant genomes. Hortic. Res. 9, uhac017 (2022).

58. Cantarel, B. L. et al. MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 18, 188-196 (2008).

59. Ray, R. et al. A persistent major mutation in canonical jasmonate signaling is embedded in an herbivory-elicited gene network. Proc. Natl. Acad. Sci. USA 120, e2308500120 (2023).

60. Li, W. & Godzik, A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658-1659 (2006).

61. Kim, D., Langmead, B. & Salzberg, S. L. HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357-360 (2015).

62. Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290-295 (2015).

63. Korf, I. Gene finding in novel genomes. BMC Bioinformatics 5, 59 (2004).

64. Lomsadze, A., Ter-Hovhannisyan, V., Chernoff, Y. O. & Borodovsky, M. Gene identification in novel eukaryotic genomes by self-training algorithm. Nucleic Acids Res. 33, 6494-6506 (2005).

65. Stanke, M., Tzvetkova, A. & Morgenstern, B. AUGUSTUS at EGASP: using EST, protein and genomic alignments for improved gene prediction in the human genome. Genome Biol. 7, S11 (2006).

66. Brůna, T., Hoff, K. J., Lomsadze, A., Stanke, M. & Borodovsky, M. BRAKER2: automatic eukaryotic genome annotation with GeneMark-EP+ and AUGUSTUS supported by a protein database. NAR Genom. Bioinform. 3, lqaa108 (2021).

67. Nevers, Y. et al. Quality assessment of gene repertoire annotations with OMArk. Nat. Biotechnol. (2024).

68. Emms, D. M. & Kelly, S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 20, 238 (2019).

69. Hu, J. et al. NextDenovo: an efficient error correction and accurate assembly tool for noisy long reads. Genome Biol. 25, 107 (2024).

70. Bi, C. et al. PMAT: an efficient plant mitogenome assembly toolkit using low-coverage HiFi sequencing data. Hortic Res. 11, uhae023 (2024).

71. Wick, R. R., Schultz, M. B., Zobel, J. & Holt, K. E. Bandage: interactive visualization of de novo genome assemblies. Bioinformatics 31, 3350-3352 (2015).

72. Tillich, M. et al. GeSeq - versatile and accurate annotation of organelle genomes. Nucleic Acids Res. 45, W6-W11 (2017).

73. Hao, Z. et al. RIdeogram: drawing SVG graphics to visualize and map genome-wide data on the idiograms. PeerJ Comput. Sci. 6, e251 (2020).

74. Nagaki, K., Kashihara, K. & Murata, M. A centromeric DNA sequence colocalized with a centromere-specific histone H3 in tobacco. Chromosoma 118, 249-257 (2009).

75. Shi, X. et al. The complete reference genome for grapevine (Vitis vinifera L.) genetics and breeding. Hortic. Res. 10, uhad061 (2023).

76. Wang, Y.H. et al. Telomere-to-telomere carrot (Daucus carota) genome assembly reveals carotenoid characteristics. Hortic. Res. 10, uhad103 (2023).

77. Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884-i890 (2018).

78. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie2. Nat. Methods 9, 357-359 (2012).

79. Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078-2079 (2009).

80. Ramirez, F. et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 44, W160-W165 (2016).

81. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, R137 (2008).

82. Vollger, M.R., Kerpedjiev, P., Phillippy, A.M. & Eichler, E.E. StainedGlass: interactive visualization of massive tandem repeat structures with identity heatmaps. Bioinformatics 38, 2049-2051 (2022).

83. Simpson, J. T. et al. Detecting DNA cytosine methylation using nanopore sequencing. Nat. Methods 14, 407-410 (2017).

84. Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094-3100 (2018).

85. Tang, H. et al. Synteny and collinearity in plant genomes. Science 320, 486-488 (2008).

86. Wang, Y. et al. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 40, e49 (2012).

87. Langdon, Q. K., Peris, D., Kyle, B. & Hittinger, C. T. sppIDer: A species identification tool to investigate hybrid genomes with high-throughput sequencing. Mol. Biol. Evol. 35, 2835-2849 (2018).

88. Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312-1313 (2014).

89. Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586-1591 (2007).

90. Nguyen, L. T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268-274 (2015).

There is NO Competing Interest.

SupplementarytablesV4.xlsx
Supplementary Tables
SupplementaryfiguresV4.docx
Extendedfigures.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

The complete genome assembly of Nicotiana benthamiana reveals genetic and epigenetic landscape of centromeres

Status:

Version 1

Abstract

Figures

Main

Results

Discussion

Declarations

References

Methods

Additional Declarations

Supplementary Files

Status:

Version 1