Gapless genome assembly and evolutionary analysis of Cnidium monnieri (Apiaceae)

Cnidium monnieri (L.) is an annual herb widely distributed across China and serves as a highly significant traditional Chinese medicine^[1]. Its fruit, known as "She Chuang Zi", has been used for centuries to treat dermatological and gynecological conditions^[1]. Modern phytochemical and pharmacological studies have revealed that the plant is rich in diverse active constituents, such as coumarins, chromones, volatile oils, glycosides, and terpenoids. A large proportion of these metabolites display notable biological activities, including antiosteoporotic, antiviral, and anti-inflammatory effects^[2−7].

Beyond its medicinal value, C. monnieri is ecologically crucial. As a nectar- and pollen-rich umbellifer, C. monnieri supports beneficial arthropods, including Coccinellidae (ladybird beetles), Chrysopidae (lacewings), and Syrphidae (hoverflies), by providing supplemental food resources and shelter^[8]. Field studies have shown that planting C. monnieri flower strips along the edges of wheat (Triticum aestivum)–maize (Zea mays) rotation fields creates "bridge habitats" that sustain ladybird beetle populations during the wheat harvest period and facilitate their migration into the adjacent maize fields, thereby enhancing aphid suppression^[9]. When infested by aphids, C. monnieri releases β-myrcene, a key herbivore-induced plant volatile (HIPV). This terpene specifically attracts the predatory ladybird beetle Harmonia axyridis by binding to the beetle's odorant-binding proteins (HaxyOBPs), thereby mediating triple-trophic plant–aphid–predator interactions and promoting natural pest regulation^[10]. These characteristics make C. monnieri an attractive candidate for ecological pest management strategies aimed at enhancing natural enemy conservation and reducing dependence on chemical pesticides.

Despite this dual importance, genomic research on C. monnieri has been limited, hindering efforts to fully elucidate the regulatory and evolutionary mechanisms underlying coumarin biosynthesis and the genetic basis of its ecological interactions. A recently published chromosome-level genome of C. monnieri^[11] provides an initial foundation; however, a gapless assembly and high-confidence genome evolution analyses are still lacking.

Here, we present a telomere-to-telomere (T2T) genome assembly for C. monnieri. By integrating complementary long-read sequencing technologies with high-throughput chromosome conformation capture (Hi-C) scaffolding, we achieved a substantial improvement in the genome's continuity and completeness. Through phylogenetic reconstruction and large-scale synteny analyses, we clarified the chromosomal evolution across the Apiaceae family. This establishes a high-quality genomic foundation for understanding trait diversification in a species of both medicinal and ecological significance.

High-quality genomic DNA was extracted and purified from the young leaves of C. monnieri via an optimized cetyltrimethylammonium bromide (CTAB) protocol. The paired-end sequencing libraries (150 bp) were prepared and subjected to sequencing on the Illumina NovaSeq X Plus platform (Illumina Inc., San Diego, CA, USA). For long-read sequencing, PacBio High-Fidelity (HiFi) and Oxford Nanopore Technologies (ONT) libraries were constructed from the young leaf tissues following the standard protocols of Pacific Biosciences Co. (Menlo Park, CA, USA) and Oxford Nanopore Technologies (Oxford, UK), respectively. The PacBio Revio platform was utilized to produce the HiFi long reads, whereas the ONT ultra-long reads were generated from the PromethION platform. To generate the Hi-C library, fresh seedling leaves of C. monnieri were harvested and treated with 2% formaldehyde under a vacuum to crosslink the chromatin. Glycine was then applied to halt the crosslinking reaction. After purification of the nuclei, chromatin was digested utilizing 100 units of the DpnII restriction enzyme, followed by biotin-14-deoxyadenosine triphosphate (dATP) labeling. The religated DNA was subsequently fragmented to a size range of 300–700 bp. These fragments underwent end-repair, A-tailing, and purification steps. Ultimately, the resulting Hi-C library was measured for concentration and sequenced on the Illumina NovaSeq X Plus platform.

Genome survey

We evaluated the genome size of C. monnieri by conducting a k-mer analysis based on the Illumina paired-end sequencing data. Specifically, KMC v3.2.1^[12] was applied to determine the frequency of the k-mers. The frequency profile (k = 17) was subsequently fed into GenomeScope v1.0^[13] to predict the overall genome size.

RNA-seq sequencing

Samples from five tissues (root, stem, leaf, flower, and fruit) were harvested from C. monnieri. Total RNA was isolated from these tissues utilizing TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). RNA sequencing libraries were then constructed according to the standard Illumina protocol and paired-end reads sequenced on the Illumina NovaSeq X Plus platform.

Genome assembly, gap-filling, and evaluation

To assemble PacBio HiFi reads, we used Hifiasm v0.19.3^[14] with the parameters '-l3 -u1'. We additionally used ONT ultra-long reads that were >100 kb long as the option '-ul'. To assemble the ONT ultra-long reads, NextDenovo v2.5.0^[15] was applied with the parameters 'read_cutoff = 10k; seed_depth = 45; nextgraph_options = -a 1 -q 16'. Chromosome-level scaffolding for both the HiFi and ONT assemblies was achieved by integrating the Hi-C data using ALLHiC v0.9.13^[16]. Finally, the interaction matrix of each chromosome was checked using HiC-Pro v3.1.0^[17]. To further fill the gap in the HiFi assembly, we mapped the ONT contigs to the pseudochromosomes of HiFi using Winnowmap v2.03^[18], and then filled the gap using the ONT contigs spanning the gap.

The final assembly quality was comprehensively evaluated. First, we ran Benchmarking Universal Single Copy Orthologs (BUSCO) v5.0.0^[19] against the embryophyta_odb10 and the viridiplantae_odb10 databases to evaluate the gene space's completeness. Then the long terminal repeat (LTR) assembly index (LAI)^[20] was computed to assess the repeat completeness using LTR_retriever v2.9.0^[21]. Finally, the assembly accuracy metrics, including the base pairs' quality value (QV) and k-mer completeness, were determined utilizing Merqury v1.3^[22].

Annotation of the transposable elements and protein-coding genes

To annotate the repeat sequences, a C. monnieri-specific repeat database was established using RepeatModeler v2.0.1^[23]. To locate repetitive elements across the genome, we ran RepeatMasker v4.0.9 (www.repeatmasker.org) against the database with the following parameters: '-e rmblast -div 40 -norna'. We used barrnap v0.9 (https://github.com/tseemann/barrnap) to identify the ribosomal DNA (rDNA) tandem repeats.

For the annotation of protein-coding genes, we used Braker v3.0.8^[24], integrating ab initio predictions with supporting evidence from both the transcriptome data and homologous proteins. The transcriptome evidence was derived from the RNA-seq data of five tissues. The raw RNA-seq reads underwent quality trimming and adapter removal via Trimmomatic v0.39^[25]. The resulting high-quality reads were mapped back to the genome with HISAT2 v2.1.0^[26], utilizing the settings '--min-intronlen 20 --max-intronlen 15,000'. The protein evidence included the protein sequences of Arabidopsis thaliana^[27], Daucus carota^[28], and the UniProt proteins of Embryophyta (www.uniprot.org). Finally, functional assignments for the predicted coding genes were performed with InterProScan v5.52-86.0^[29].

Phylogenetic analysis and estimation of divergence times

To resolve the evolutionary relationships within the Apiaceae family, we used protein sequences from C. monnieri and 10 previously reported Apiaceae species, including Centella asiatica^[30], D. carota^[28], Cryptotaenia japonica^[31], Oenanthe sinensis^[32], Angelica sinensis^[33], Apium graveolens^[34], Heracleum sosnowskyi^[35], Coriandrum sativum^[36], Peucedanum praeruptorum^[37], and Saposhnikovia divaricata^[38] (Supplementary Table S1). The outgroups were A. thaliana^[27] and Vitis vinifera^[39]. OrthoFinder v2.5.4^[40] was applied to cluster the orthologous gene families. Following this, the protein sequences of identified single-copy genes from each species were extracted. These sequences were aligned using MAFFT v7.487^[41] and subsequently trimmed by removing poorly aligned regions with trimAl v1.5.0^[42]. The aligned sequences were then used to reconstruct a maximum likelihood phylogenetic tree via IQ-TREE v3.0.1^[43], accompanied by an ultrafast bootstrap analysis consisting of 1,000 replicates (-bb 1,000). To infer the temporal divergence of these lineages, r8s v1.81^[44] was utilized, incorporating fossil calibration points retrieved from the TimeTree database^[45].

Comparative synteny analysis

For evaluating the chromosomal synteny, collinearity among different species was examined via mcscan (JCVI v1.1.19^[46]), utilizing orthologous single-copy genes derived from OrthoFinder. The syntenic blocks between species were identified by applying MCScanX^[47], setting a threshold of at least 10 matching genes per syntenic block (-s 10 -m 50 -b 2). The resulting synteny relationships were subsequently plotted with the SynVisio (https://synvisio.github.io)^[48].

The analysis of the karyotype indicated that C. monnieri was a diploid species containing 20 chromosomes (2n = 2x = 20) (Supplementary Fig. S1). According to the k-mer analysis, the approximate size of the C. monnieri genome was 1.17 Gb (Supplementary Fig. S2). To construct a high-quality genome assembly for C. monnieri, we performed PacBio HiFi long reads and ONT ultra-long reads sequencing. We generated ~99.7 Gb PacBio HiFi long reads (read N50 = 18.3 kb) and ~111.2 Gb ONT ultra-long reads (read N50 = 100.9 kb), which covered ~85× and ~95× of the C. monnieri genome, respectively (Supplementary Tables S2 and S3). Initial assembly of the PacBio HiFi data was conducted using Hifiasm^[14], yielding 512 contigs with an N50 of 100.9 Mb (Supplementary Table S2). We also assembled the ONT ultra-long reads using NextDenovo^[15], resulting in 670 contigs with an N50 of 5.8 Mb (Supplementary Table S2). Since the contigs assembled from the HiFi reads exhibited a notably higher N50 than those assembled from the ONT reads, we adopted the HiFi contigs for chromosome scaffolding (Supplementary Fig. S3). These contigs were anchored, oriented, and manually curated into 10 pseudochromosomes using Hi-C data with ALL-HiC^[16]. The resulting assembly contained only one gap, and the anchoring rate of the HiFi contigs was 95.9% (Supplementary Table S2), indicating the high continuity of the HiFi assembly. We then filled the gap by using the ONT contigs spanning the gap (Supplementary Fig. S4). Finally, we performed telomere patching using TeloClip v0.3.2 (https://github.com/Adamtaranto/teloclip) along with raw long reads, yielding all 20 telomeres that were enriched with seven-base telomere repeats (5'-TTTAGGG-3') (Supplementary Fig. S5; Supplementary Table S4). The total size of the pseudochromosomes in the final assembly was 1.13 Gb (Fig. 1a, b and Supplementary Table S2), which is consistent with the estimated genome size or that of Cmo_YZ (1.14 Gb)^[11]. Collinearity analysis with the Cmo_YZ assembly revealed that the current assembly showed high collinearity with it, though there are several chromosomal inversions (Fig. 1c).

Figure 1. 

Assembly of the C. monnieri genome. (a) The circos plot of features across the chromosomes of C. monnieri. Tracks i–vi represent the chromosomes, LTR/Gypsy, LTR/Copia, DNA-TEs, gene, and GC content, respectively. (b) Heatmap displaying whole-genome Hi-C interaction frequencies for C. monnieri. The intensity of the interactions is represented by colors shading from yellow (low) to red (high). (c) Collinearity mapping between this study (top) and Cmo_YZ^[11] (bottom). Gray ribbons connect collinear regions between the two genomes. Gap regions are represented by yellow blocks in Cmo_YZ assembly. The triangles at the ends of chromosomes indicate the identified telomeric repeat sequences.

To characterize the centromeres, we performed a comprehensive identification of tandem repeats across the genome using TRF^[49] and identified two predominant repeat units (118 and 180 bp of monomer lengths, respectively; Supplementary Fig. S6). The 118-bp unit (Cmon118) was the potentially primary centromeric satellite repeat. Although this repeat dominated the centromeric regions of seven chromosomes, it was absent on chromosomes 5, 9, and 10 (Supplementary Fig. S7). This may be because the centromeres of these three chromosomes are maintained through the massive accumulation of LTR/Gypsy retrotransposons^[50]. We also discovered another highly abundant 180-bp repetitive sequence (Cmon180). Compared with Cmon118, this unit was highly conserved and was predominantly distributed in the terminal regions of the chromosomes, indicating that it was likely a subtelomeric satellite repeat^[51].

To assess the quality of the assembly, we first evaluated its completeness using BUSCO analysis with embryophyta_odb10 and viridiplantae_odb10, of which 98.2% (1,584/1,614) and 98.8% (420/425) complete genes were present in the assembly, respectively (Supplementary Table S5). Additionally, we estimated the contiguity and repeat region assembly quality using the LAI^[20], resulting in a value of 20.94 (Supplementary Table S2), indicating that the current assembly had high completeness. We also mapped the PacBio HiFi reads, Illumina reads, and RNA-seq reads to the assembly, yielding mean mapping rates of 99.94%, 99.78%, and 94.27%, respectively. Further evaluation using Merqury^[22] showed that the average QV score of the chromosomes was 55.86 (ranging from 54.94 to 57.29 for each chromosome), indicating high (> 99.999%) assembly accuracy (Supplementary Table S6). Additionally, the Hi-C interaction matrices also showed that the assembly had high consistency (Fig. 1b and Supplementary Fig. S7). Compared with the previously published C. monnieri genome (Cmo_YZ^[11]), our assembly achieved a T2T status. Consequently, our assembly exhibits a higher contig N50 (100.9 Mb vs. 78.45 Mb) and chromosome anchoring rate (95.9% vs. 93.9%). Overall, our gapless assembly represents a contiguous and complete genome of C. monnieri.

Annotation of TEs and protein-coding genes in the C. monnieri genome

In total, 884.86 Mb of repeat elements was annotated, constituting 74.76% of the C. monnieri genome (Supplementary Table S7). Among the annotated transposable elements (TEs), LTR retrotransposons were the most dominant TE type, covering 51.45% of the genome sequence (Fig. 2a, Supplementary Table S7). In contrast, the content of other TEs (e.g., DNA transposons, long interspersed nuclear elements (LINEs), short interspersed nuclear elements (SINEs), and rolling-circle (RC)/Helitron) was low (< 3%) (Fig. 2a and Supplementary Table S7). Analysis of the distribution of sequence divergence across the TE superfamilies revealed that younger TEs (with lower divergence), particularly for LTR retrotransposons, had higher proportions, indicating the LTR retrotransposons had recently been active in the genome (Fig. 2a). Notably, we found a strong positive correlation between genome size and TE content in Apiaceae species (R = 0.98, p = 3.35 × 10⁻⁶; Supplementary Fig. S8).

Figure 2. 

Characterization of TEs and protein-coding genes in the C. monnieri genome. (a) A repeat landscape of the C. monnieri genome. The x-axis (divergence) represents the Kimura 2-parameter distance between individual repeat elements and their corresponding consensus sequences. The inset pie chart illustrates the percentage of TEs within the whole genome. (b) BUSCO completeness for the different datasets. (c) Venn diagram of the annotation results from different databases.

Subsequently, we predicted the protein-coding genes from the repeat-masked genome using Braker^[24], which integrates ab initio gene predictions, transcriptome data, and homologous protein evidence. In total, 39,218 protein-coding genes were annotated in the C. monnieri genome, of which 99.06% were located on the chromosomes. We then evaluated the completeness of gene annotation using BUSCO^[19], resulting in 99.6% (1,607/1,614 in embryophyta_odb10) and 99.7% (424/425 in viridiplantae_odb10) of all BUSCO genes being contained in the annotated gene set (Fig. 2b, Supplementary Table S8). Of all the annotated genes, 97.48% (38,231/39,218) had functional annotations from the InterPro database^[29], and 87.79%, 71.74%, 58.65%, and 56.06% contained functional terms in the databases of PANTHER, Pfam, Gene3D, and SUPERFAMILY, respectively (Fig. 2c, Supplementary Table S9).

Genome evolution of Apiaceae

To resolve the evolutionary history of Apiaceae and elucidate the chromosomal diversification patterns within the clade, we reconstructed a phylogenetic tree using 321 single-copy orthologous genes across the Apiaceae species with A. thaliana and V. vinifera as the outgroups (Fig. 3a). The resulting maximum likelihood tree was used to infer the divergence times via fossil-calibrated molecular clock methods. Consistent with previous studies, species belonging to the subfamily Apioideae formed a well-supported monophyletic group, and C. monnieri was resolved as the sister to P. praeruptorum and S. divaricata^[37]. C. monnieri, along with P. praeruptorum and S. divaricata, diverged from C. sativum 11 million years ago (Mya) (Fig. 3a).

Figure 3. 

Genome evolution of C. monnieri and Apiaceae. (a) Phylogenetic analysis of C. monnieri and other Apiaceae lineages. (b) Syntenic analysis of the Apiaceae species. (c) Synteny blocks between A. sinensis and P. praeruptorum. (d) Synteny blocks between A. sinensis and C. monnieri. For (c) and (d), the colored ribbons represent syntenic blocks identified by MCScanX and visualized using SynVisio. Numbers indicate chromosome identification numbers (IDs).

To investigate the patterns of chromosomal evolution, we conducted a genome-wide synteny analysis based on single-copy ortholog anchors across Apiaceae species, revealing the uneven conservation of structure of chromosomes in Apiaceae (Fig. 3b). Species with 11-chromosome karyotypes in Apioideae, including A. sinensis, H. sosnowskyi, A. graveolens, C. sativum, and P. praeruptorum, exhibited highly conserved syntenic blocks (Fig. 3b). The detailed synteny analyses between A. sinensis (as the reference) and other 11-chromosome species further resolved fine-scale chromosomal dynamics (Supplementary Fig. S9). A. sinensis consistently displayed a low degree of chromosomal structural rearrangement across the comparisons, indicating it preserved one of ancestral and conserved karyotypes within the Apioideae. For example, the genome of A. graveolens showed only limited rearrangement relative to A. sinensis, with chromosomal rearrangements detected only for chromosomes 9, 10, and 11 (Supplementary Fig. S9a). As the phylogenetic distance increased, the number of structural rearrangements generally rose; however, even in more distantly related 11-chromosome species, certain chromosomes remained highly conserved, such as chromosome 2 of P. praeruptorum (Fig. 3c), representing conserved linkage blocks that had persisted throughout diversification of the Apioideae. This is consistent with prior studies indicating that the ancestral chromosome number of Apioideae was 11^[34,52], suggesting that these species retain the plesiomorphic chromosomal architecture of the subfamily. In contrast, species with deviant chromosome counts exhibited markedly reduced collinearity and increased interchromosomal rearrangements (Supplementary Fig. S10). For instance, C. monnieri (n = 10) underwent extensive chromosomal rearrangements relative to 11-chromosome species like A. sinensis, and all chromosomes except 7 and 10 experienced extensive chromosomal rearrangements, leading to its reduced chromosome count (Fig. 3d).

This pattern was further reflected in the comparative quantification of interchromosomal rearrangements by measuring how many chromosomal segments from other species contributed to a target chromosome. Species with 11 chromosomes (e.g., A. sinensis and A. graveolens) had fewer such segments, indicating lower rearrangement frequencies (Supplementary Fig. S11). A notable feature was the genome of O. sinensis (n = 21), which underwent an additional whole-genome duplication (WGD)^[32]. A detailed chromosome-to-chromosome synteny comparison showed that most chromosomes of O. sinensis matched two counterparts in A. sinensis, providing strong evidence for a recent WGD event. Nevertheless, one set of its chromosomes still maintained the 11-chromosome structure, which explains why it had the fewest interchromosomal rearrangements (Supplementary Fig. S10d).

These results demonstrated that WGD, conservation of the ancestral 11-chromosome architecture, and lineage-specific rearrangements collectively drove the genomic evolution of Apiaceae. Species with 11 chromosomes in the Apioideae represent a conserved lineage, while C. monnieri and other taxa with deviant chromosome numbers have undergone extensive rearrangements. The observed genomic dynamics shed light on the adaptive evolution and taxonomic diversification of the Apiaceae.

In this study, we generated a gapless T2T genome assembly of C. monnieri, achieving a significant improvement in continuity and accuracy over previously reported assemblies^[11]. The high QV score, nearly complete BUSCO coverage, and the identification of 20 telomeres highlighted the high quality and structural completeness of the genome. This T2T genome provides a robust foundation for downstream research on medicinal compound biosynthesis, functional genomics, and evolutionary biology. However, limitations inherent to the haplotype-collapsed assembly strategy remain. For example, chromosome 7 displayed the signatures of potential chimeric assembly, a phenomenon most likely caused by the extreme sequence divergence between the two haplotypes of this chromosome. Generating a completely phased diploid T2T genome to definitively resolve these complex allelic variations will be the primary focus of future research.

The repeat composition of the C. monnieri genome revealed a substantial accumulation of LTR retrotransposons, which is consistent with patterns reported in several Apiaceae species^[36−38]. The relatively high proportion of LTRs, particularly Gypsy and Copia elements, suggests that TE expansion has played a major role in shaping variations in genome size within the family, supporting the hypothesis that differences in LTR proliferation and deletion rates may underlie the pronounced genome size disparities observed among closely related genera^[53].

Comparative genomic analyses further clarified the evolutionary placement of C. monnieri within the Apiaceae. Using 321 conserved single-copy orthologs, we reconstructed a robust phylogeny that positioned C. monnieri as sister to P. praeruptorum and S. divaricata. The estimated divergence times among these species provide new insights into lineage diversification patterns within the tribe and highlight the value of high-quality genomic datasets for reconstructing Apiaceae's evolution. Our comparative genomic analyses illuminated substantial differences in chromosomal evolution across the Apiaceae. Species retaining the ancestral 11-chromosome karyotype, such as A. sinensis and A. graveolens^[33,34], showed high syntenic conservation, whereas C. monnieri (n = 10) displayed extensive chromosomal rearrangements involving both fissions and fusions. These lineage-specific structural changes may underlie unique metabolic traits, ecological functions, and adaptation strategies.

Overall, our study offers a significantly improved genomic resource for C. monnieri, facilitating future multidisciplinary research. First, from a pharmacological perspective, C. monnieri is rich in active secondary metabolites, particularly coumarins and volatile oils. High-continuity genomic data are crucial for the precise identification of biosynthetic gene clusters and provide a reliable foundation for elucidating the complete biosynthetic pathways of these medicinal compounds. Second, from an ecological standpoint, C. monnieri mediates triple-trophic plant–aphid–predator interactions by releasing specific HIPVs such as β-myrcene. The complete genome sequence will greatly facilitate the discovery of the genetic mechanisms underlying these ecological traits and volatile emissions. Finally, integrating this T2T assembly with future pan-genomic analyses and transcriptomic profiling will serve as a robust resource for population genomics and marker-assisted breeding. This will ultimately accelerate the cultivation of superior C. monnieri varieties with enhanced medicinal yields and ecological pest management capabilities, further advancing both the biological and applied value of this species.