C. longepaniculata leaves were obtained in the Chinese province of Sichuan (27°50’ N, 105°20’ E). Genomic DNA was then extracted and processed according to the recommended library preparation methodology by PacBio. Subsequently, the DNA was utilized to construct PacBio sequencing libraries using the standard single molecular real-time bell construction protocol provided by PacBio. The libraries were then sequenced on the PacBio CCS platform. Adapter sequences and short, low-quality reads were removed from the PacBio HiFi reads using SMRTLink (parameters: –min-passes = 3 –min-rq = 0.99). Using the Agencourt AMPure XP-Medium Kit, the same DNA was used to create an Illumina paired-end library with insert sizes ranging from 200 to 400 bp. The Illumina NovaSeq 6000 platform was then used to sequence the resulting library. On the MGISEQ2000 platform, a Hi-C library was sequenced in accordance with the proximo Hi-C plant methodology (https://info.phasegenomics.com/protocols) developed by Phase Genomics (Seattle, WA, United States). After removing adaptors and duplicate reads, raw reads were filtered with a sequencing quality of >Q30 ( Supplementary Table S1 ).

Lastly, RNA sequencing was performed on a leaf sample obtained from the same individual. Using the Qiagen RNeasy Plant Mini kit, total RNA was extracted from each sample. Then, with the TruSeq RNA Library Preparation kit, RNA-seq libraries were created, and 11.74 Gb of RNA-seq sequences were obtained by sequencing on the Illumina NovaSeq 6000 platform (Belton et al., 2012) ( Supplementary Table S1 ).

Utilizing FASTP v0.20.0 (Chen et al., 2018), the short-read sequences were filtered to exclude adaptor contamination and low-quality reads. Subsequently, genome size, repeat content, and heterozygosity rate were estimated using GCE v1.0.2 tool from 19-mer histograms (Liu et al., 2013).

With the parameters “‘–write-ec –write-paf -u -l 0’’, the PacBio HiFi reads were assembled using Hifiasm v.0.15 (Cheng et al., 2021) and redundant sequences were removed using Purge Haplotigs (Roach et al., 2018). Hi-C data was used to anchor the contigs to the chromosomes. In summary, genuine interaction pairs were obtained by filtering the Hi-C reads using FASTP v0.20.0 ( Supplementary Table S5 ; Chen et al., 2018) and selecting unique mapped read pairs using the HiCUP v0.8.0 (Wingett et al., 2015) pipeline. Next, a chromosome-level assembly was created using 3D-DNA v201013 (https://github.com/aidenlab/3d-dna), after the sequence had been aligned against the draft genome assembly using Juicer (Dudchenko et al., 2017). The “embryophyta_odb10” ortholog set was used to assess the quality of the genome using BUSCO v4.05 (Simão et al., 2015).

To ensure the integrity of the C. longepaniculata genome assembly, high-quality short-read sequences were mapped back to the genome using the BWA v0.7.17 (Li and Durbin, 2010). This process helps to verify the correct alignment and orientation of the assembled sequences. Following the mapping, duplicate reads were removed using Picard v2.23.6 (broadinstitute.github.io/picard) to prevent redundancy that could skew variant identification results. Subsequently, SNPs (Single Nucleotide Polymorphisms) and InDels (insertions and deletions) were identified using SAMtools v1.10 (Li et al., 2009), which provides insights into potential genomic variations and polymorphisms.

In addition to the above steps, HiFi reads were mapped back to the genome using Minimap2 (Li, 2018). This step further ensures that the high-fidelity, long-read sequencing data aligns accurately with the assembled genome, confirming the assembly’s overall quality and completeness.

Repetitive sequences in the de novo repeat library were predicted using RepeatMasker v4.0.7 (Tarailo-Graovac and Chen, 2009), and a de novo repeat library was created using RepeatModeler (Flynn et al., 2020) and LTRfinder v1.07 (Xu and Wang, 2007) with default settings. In parallel, homologous repeat prediction was carried out using the Repbase v21.12 database and RepeatMasker v4.0.7 and RepeatProteinMasker v4.0.7 (http://www.repeatmasker.org/cgi-bin/RepeatProteinMaskRequest) (Bao et al., 2015). To create nonredundant repeated sequences, the two sets of anticipated repetitions were concatenated. With Tandem Repeat Finder v4.09 (Benson, 1999), tandem repetitions were found.

Gene prediction in the repeat-masked genome was done using three different methods: homology searching, reference-guided transcriptome assembly, and ab initio prediction. According to Keilwagen et al. (2016), homologous peptides from related species (Cinnamomum kanehirae, Phoebe bournei, Litsea cubeba, Phoebe kanehirae, and Phoebe zhennan) were aligned with the assembly using GeMoMa v1.6.1 (Keilwagen et al., 2016). Gene structure information was then collected for homolog prediction. Filtered mRNA-seq data were used to align STAR v2.7.3a (Dobin et al., 2013) to the reference genome in order to perform RNA-seq-based gene prediction (Dobin et al., 2013). PASA v2.3.3 (Haas et al., 2008) was used to predict open reading frames (ORFs), and StringTie v1.3.4d (Pertea et al., 2015) was used to assemble the transcripts. A training set was also produced for the de novo prediction. The training data was then used for ab initio gene prediction using Augustus v3.3.1 (Stanke et al., 2008) with default settings. After creating an integrated gene collection with EvidenceModeler (EVM, v1.1.1; Haas et al., 2008), miscoded genes were further filtered and genes containing TEs were eliminated using the TransposonPSI program (http://transposonpsi.sourceforge.net/) (Urasaki et al., 2017). PASA was used to identify alternative splicing regions and untranslated regions (UTRs) based on RNA-seq assemblies. For each locus, we kept the longest transcripts, and areas outside of the ORFs were called UTRs. By examining five protein/function databases, the projected protein-coding genes’ domains, motifs, and gene function information were found. Predicted protein-coding genes were thoroughly annotated using InterProScan v5.36 (Zdobnov and Apweiler, 2001), which included predicted signal peptides, transmembrane topologies, functional classifications, GO keywords, protein motifs and domains, and protein family identification. KO words were found using KaaS (https://www.genome.jp/kegg/kaas/) searching the KEGG database. A 1e-5 E value limit was utilized for searches against the Swiss-Prot (Bairoch and Apweiler, 2000), NR (Marchler-Bauer et al., 2011), TrEMBL (Bairoch and Apweiler, 2000), and COG databases (Galperin et al., 2015) using BLASTP v.2.7.1 (Altschul et al., 1997). From these database searches, the top hits were combined to get the final results.

Database searches and model-based prediction were the two methods employed to find noncoding RNA sequences in the genome. With eukaryote parameters, tRNAscan-SE v2.0.12 (Lowe and Eddy, 1997) was used to predict tRNAs. INFERNAL v1.1.5 (Nawrocki and Eddy, 2013) was utilized to search the Rfam database and identify sequences of microRNA, rRNA, small nuclear RNA, and small nucleolar RNA. RNAmmer v1.2 (Lagesen et al., 2007) was used to forecast the rRNAs and their components.

In order to investigate the evolutionary background of C. longepaniculata, we obtained the genomes of eleven different species from publicly accessible databases. These species include Aquilegia coerulea, Arabidopsis thaliana, Amborella trichopoda, Liriodendron chinense, Nymphaea colorata, Piper nigrum, Cinnamomum kanehirae, Persea americana, Chimonanthus salicifolius, and Litsea cubeba ( Supplementary Table S2 ). Furthermore, genes with glaring mistakes were eliminated and the longest transcripts were selected as samples for genes having alternative splice variants. The homology matrix of the orthogroups (gene families) among these chosen species was also inferred using OrthoFinder v2.5.4 (Emms and Kelly, 2019). For the purpose of creating phylogenetic trees using the maximum-likelihood (ML) approach, 230 single-copy gene groupings were found. The single-copy gene groups were subjected to multiple sequence alignment using MAFFT v.7.453 (Katoh and Standley, 2013). Subsequently, the coding sequences (CDS sequences) were aligned using PAL2NAL v.14 (Suyama et al., 2006) in accordance with the alignments of the respective proteins. Each single-copy gene group was given a maximum likelihood phylogenetic tree using IQ-TREE v.2.1.4-beta (Nguyen et al., 2015). Ultimately, these 230 single-copy gene trees were combined using the ASTRAL v.5.6.3 (Zhang et al., 2022) software to create a species tree based on the multispecies coalescent model. We employed a Bayesian relaxed molecular clock technique to estimate the divergence periods between species using the MCMCTree tool in the PAML v.4.9 (Yang, 2007) package. Furthermore, we acquired four fossil constraints for divergence time calibration from the TimeTree website (http://www.timetree.org): Between A. trichopoda and N. colorata, 179–205 Mya; between A. coerulea and P. nigrum, 151–170 Mya; between A. coerulea and A. thaliana, 126–133 Mya; and between A. thaliana and P. tremuloides, 102–113 Mya. Moreover, gene families in phylogenetic tree-constructing species that had experienced expansion or contraction were identified using CAFÉ v.4.2.1 (De Bie et al., 2006). Genes belonging to particular extended gene families were then functionally analyzed using GO and KEGG enrichment with TBtools (Chen et al., 2023).

To find members of the gene families implicated in the pathways leading to terpen biosynthesis, we conducted searches using BLASTP and Hidden Markov Models (HMMs). In particular, we gathered genes that are known to be connected to this A. thaliana pathway. Then, BLASTP (e < 1e-5) was performed to discover pathway genes in the genome of C. longepaniculata using these genes as query sequences. HMM data for each gene family’s conserved protein domain were downloaded in the interim using the Pfam website (https://pfam.xfam.org/). These HMM files were used in batch searches, along with HMMER v3.2.1 (Johnson et al., 2010). These HMM files served as the basis for batch searches. The BLASTP-identified candidate genes that lacked the matching domain were eliminated. We utilized the PF02458 Pfam domain for BAHD searching and the PF01397 and PF03936 Pfam domains for HMMER searching in order to identify the TPS genes. In order to differentiate between distinct subfamilies, phylogenetic trees were created for each gene family’s candidate genes using IQTREE v.2.1.4-beta (Nguyen et al., 2015).

We used transcriptome data from the leaf to do RNA-seq analysis to look into the expression levels of genes involved in the terpene synthesis pathway. We first performed a raw data quality assessment for transcriptome analysis using FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Subsequently, quality control was performed using FASTP v0.20.0 to filter out raw reads containing poly-N sequences and low-quality reads (< Q30). We used the high-quality reads obtained from the quality control step for the subsequent expression analysis. After filtering the high-quality RNA-seq reads, we used HISAT2 v.2.1.0 (Kim et al., 2015) to map them to the genome of C. longepaniculata. Finally, using StringTie v.1.3.4d (Pertea et al., 2015), we quantified gene expression levels in transcripts per kilobase million (TPM).

In order to create the genome of C. longepaniculata ( Figure 1A ), we initially generated 49.43Gb of paired-end reads using the Illumina NovaSeq 6000 ( Supplementary Table S1 ). Our 19-mer analysis indicated an estimated genome size of 702 Mb, with a heterozygosity rate of 3.31% ( Supplementary Table S3 ). Building on this, we utilized 31.54 Gb of high-quality high-fidelity (HiFi) reads produced on the PacBio platform, achieving 44.91-fold coverage of the C. longepaniculata genome ( Supplementary Table S1 ). Utilizing the Hifiasm tool, this sequencing effort yielded thirty contigs, totalling 768.10 Mb with an N50 size of 41.12 Mb ( Supplementary Table S4 ). To enhance the assembly, we incorporated Hi-C data, employing Hi-C-Pro v2.8.1 to recognize and retain 89.18 Gb of valid interaction paired reads ( Supplementary Table S5 ). This integration allowed the assembly to be further consolidated. Ultimately, we anchored 766.69 Mb (99.82%) of the contig sequences into 12 chromosomes ( Figures 1B, C ; Supplementary Table S6 ). The final scaffold N50 was improved to 63.78 Mb, with the longest scaffold measuring 111.09 Mb ( Table 1 ; Supplementary Table S6 ). Additionally, a heatmap of chromosome interactions demonstrated the completeness and robustness of the genome assembly ( Figure 1C ).

A combination of homology-based, reference-guided transcriptome assembly, and ab initio gene approaches were used to predict gene models in the C. longepaniculata assembly. The consensus gene set was then generated by combining the gene prediction results with EVM software. To improve the quality of gene prediction, genes with transposable elements and miscoded genes were eliminated. In the end, 39,173 genes made up the final gene list that we acquired. These protein-coding genes had an average gene length of 7,152.40 bp and an average CDS length of 1,146.30 bp ( Supplementary Tables S7, S8 ). There were 4.38 exons on average per gene, with an average exon length of 261.80 bp and intron length of 1608.33 bp ( Supplementary Table S8 ).

Approximately sixty percent of the whole-genome assembly consisted of the detected repetitive sequences (463.83 Mb) ( Supplementary Table S9 ). The most common repeat types were DNA elements (19.72%) and long terminal repeat (LTR) retrotransposons (31.96%) ( Supplementary Table S9 ). LINE (3.19 Mb) made up 0.42% of the whole genome assembly, in contrast. With computed average lengths of 106.89, 74.99, 191.27, and 118.18 bp, respectively, 6,531 miRNAs, 582 transfer RNAs (tRNAs), 3,312 rRNAs, and 634 snRNAs were found ( Supplementary Table S10 ).

With short reads an overall mapping rate of 99.20% was obtained, encompassing 99.92% of the assembly ( Supplementary Table S11 ). Moreover, HiFi reads were mapped back using Minimap2 (Li, 2018), resulting in an assembly coverage of 100% and an overall mapping rate of 99.92% ( Supplementary Table S11 ). These findings imply that the genetic information provided in our assembly was almost entirely complete. Moreover, SAMtools v1.4 (Li et al., 2009) was used to identify and filter single-nucleotide polymorphisms (SNPs), resulting in the identification of 11,919,258 heterozygous SNPs, 18,760 homozygous SNPs, and 9,182 homozygous INDELs with 5× sequencing depth. The high precision of the assembly (99.996385%) is supported by the low frequencies of homozygous SNPs and INDELs, which together make up 0.002442% and 0.001195% of the assembled genome, respectively. Ultimately, a scatter plot of the sequencing depth vs the GC-content based on 10-kb windows showed that the completed C. longepaniculata genome was free of foreign DNA contamination ( Figure 1B , Supplementary Figure S1 ). Additionally, we used the embryophyta odb10 database (https://busco.ezlab.org/) to perform BUSCO analysis on the data (Simão et al., 2015). Out of 1,614 conserved plant genes, 4.52% of them contained duplicates, yet 99.01% of them had full coverage in the genome, 0.43% were fragmented, and only 0.56% were missing ( Supplementary Table S12 ). These results clearly show that the genome assembly of C. longepaniculata we have produced is of excellent quality and has the potential to be useful for future research. The completeness of these anticipated genes was further evaluated using BUSCO analysis, yielding a BUSCO score of 93.31% (single = 89.34%, duplicated = 3.97%, fragmented = 4.96%, missing = 1.73%, genes = 1,614; Supplementary Table S13 ). Furthermore, functions for 65.80% (25,777), 35.86% (14,049), 71.29% (27,928), 36.48% (14,290), 92.78% (36,345), and 93.26% (36,533) of the genes were found by searching the Swiss-Prot, KEGG, COG, GO, TrEMBL, and NR databases, respectively ( Supplementary Table S14 ). A total of 36,632 protein-coding genes, or 93.51% of them, had their conserved functional motifs or functional terms effectively annotated ( Supplementary Table S14 ). These findings suggest that,the C. longepaniculata genome’s annotated gene set is rather comprehensive.

As described in the Materials and Methods section, 288 single-copy orthologous genes were discovered by analyzing the genomes of 12 different species. C. longepaniculata and C. camphora were found to have a closer relationship based on phylogenetic analyses of these genes. C. longepaniculata split from C. camphora and C. kanehirae between 6.01 million and 62.19 million years ago, according to projected divergence times ( Figure 2A ). Furthermore, the collinearity findings revealed a high level of collinearity among the three species ( Figure 2B ), and C. longepaniculata was discovered to have a WGD event shared by Lauraceae species. In addition, C. longepaniculata shared 10,384 gene families with four other species (C. camphora, C. kanehirae, C. salicifolius, and P. americana) and had 1,687 unique gene families ( Figure 2C ). We further investigated structural variations (SVs) among C. longepaniculata, C. camphra, and C. kanehhirae. We annotated a total of 695,731 SVs (including small insertions and deletions) between C. longepaniculata and C. camphra, and 801,107 SVs between C. longepaniculata and C. kanehhirae after whole-genome alignment using minimap2 and SV calling using SyRI ( Supplementary Figure S3 ). After filtering, we focused on SVs longer than 10kb, resulting in the retention of 949 (between C. longepaniculata and C. camphora) and 701 (between C. longepaniculata and C. kanehirae) SVs, respectively. The analysis of functional enrichment of genes within 10kb of SV breakpoints showed significant enrichment in terms of terpene synthase activity (GO:0010333), terpenoid biosynthetic process (GO:0016114), monoterpene metabolic process (GO:0016098), and chloroplast RNA processing (GO:0031425), among others. This suggests that post-WGD diploidization may play a part in improving the terpenoid biosynthesis pathway ( Supplementary Table S15 ).

After a thorough evaluation of gene families in each of the 12 species, it was determined that 1,769 extended families and 1,884 contracted families are included in the genome of C. longepaniculata ( Figure 2A ). Statistical analysis revealed that 584 expanded and 497 contracted gene families were statistically significant (P < 0.05; Supplementary Tables S16 , S17 ). These gene families were clarified by the Gene Ontology (GO) enrichment study. Significantly, gene families linked to the following were found to have expanded: “transferase activity (GO: 0016740)”, “catalytic activity (GO: 0003824)”, “response to oomycetes (GO: 0002239)”, “regulation of hydrogen peroxide metabolic process (GO: 0010310)”, and “secondary metabolic process (GO: 0019748)” ( Supplementary Table S16 ). These gene families may help to explain the high levels of oil content and environmental resistance in C. longepaniculata. The gene families for “vesicle coating (GO: 0006901), Golgi vesicle transport (GO: 0048193), response to hydroperoxide (GO: 0033194), and channel activity (GO: 0015267)” on the other hand, showed shrinkage ( Supplementary Table S17 ). Moreover, processes like “amino acid activation (GO: 0043038),” “peptide metabolic process (GO: 0006518),” “calcium-mediated signaling (GO: 0019722),” “organic acid metabolic process (GO: 0006082),” and “small molecule metabolic process (GO: 0044281)” were identified by GO enrichment analysis focused on C. longepaniculata specific gene families. These processes are all crucial to the process of oil biosynthesis ( Figure 2C ; Supplementary Table S18 ).

Due to the high concentration of terpenoids (volatile organic compounds) in the leaf essential oil of C. longepaniculata and their role as important sources of plant fragrance, we focused on the terpene biosynthesis pathway genes in the C. longepaniculata genome. Terpenoids (monoterpenes, sesquiterpenes and iridoids) are usually synthesized via the MVA and MEP pathways ( Figure 3A ). We identified a total of 37 relevant genes from these two pathways in the C. longepaniculata genome ( Supplementary Table S19 ). The results showed that there was no significant difference in the copy number of these genes in the Lauraceae species overall. However, the MVD genes, which encodes diphosphomevalonate decarboxylase, showed significant expansion compared to other species, reaching up to 5 copies in C. longepaniculata ( Supplementary Table S19 ). These genes are responsible for catalyzing the conversion of mevalonate-5-diphosphate to isopentenyl diphosphate (IPP), which is one of the important precursors for synthesizing various terpenoids. This reaction step links the early steps of terpene biosynthesis with the subsequent conversion steps. The expansion of the MVP gene family in the genome of C. longepaniculata may be of significance for its terpenoid biosynthesis.

On the other hand, 1,8-Cineole, α-terpineol, and γ-terpinene are important constituents of the leaf essential oil in C. longepaniculata, and these three compounds all belong to the class of monoterpenes. Subsequently, we focused on terpene synthase (TPS), which is crucial in the biosynthesis of terpenes. It catalyzes a variety reactions in the MVA and MEP pathways to produce the basic skeletons of terpenoids. We identified 74 genes encoding terpene synthase in C. longepaniculata ( Supplementary Table S20 ). Notably, the TPS gene family in Lauraceae species significantly expanded compared to other selected species (C. longepaniculata: 74; C. camphora: 75; C. kanehirae: 88; Supplementary Table S20 ). More tree-based study of these TPS genes showed that the TPS-b and TPS-g subfamilies have grown in Lauraceae species ( Figure 3B ; Supplementary Table S20 ). Members of these two subfamilies mainly synthesize non-cyclic monoterpenes, sesquiterpenes, and diterpenes. This partially clarifies the significant contribution of the TPS-b and TPS-g gene families’ expansion to the accumulation of monoterpenes in C. longepaniculata’s leaf essential oil. BAHD acyltransferases, which participate in the synthesis of various flavors and fragrances in plants, generally transform terpenes into esters. In the genome of C. longepaniculata, we identified 98 genes that belong to the BAHD gene family. Phylogenetic analysis showed that the growth of different subfamilies of the BAHD gene family was significantly different between species ( Figure 3C ; Supplementary Table S21 ). We observed expansion trend in several subfamilies of Lauraceae species, including BAHD-I-b, BAHD-II-b, and BAHD-III-a. Specifically, for C. longepaniculata, the number of genes in the BAHD-II-b subfamily reached up to 19 ( Supplementary Table S21 ). Also, the BAHD-II-b and BAHD-III-a subfamilies of C. longepaniculata formed a separate sub-branch compared to the reference sequence of A. thaliana. This suggests that the lineage specificity of the BAHD-II-b and BAHD-III-a genes may have different functions, which may help create the unique smell of C. longepaniculata. Observation of the chromosomal locations of TPS and BAHD gene families revealed that, similar to previous results in other plants (Xu et al., 2022), many TPS genes in C. longepaniculata exhibited tandem repeats. Specifically, these TPS genes formed gene clusters on chromosomes 3, 4, 7, and 9, showing a distinct clustered distribution pattern. On chromosomes 7 and 9, we found two large gene clusters, one with 18 genes and the other with 27 genes. This suggests that tandem duplication events recently created these genes ( Figure 3D ). At the same time, we observed several gene clusters on chromosomes 2 and 10 for the distribution of BAHD genes their overall distribution did not follow the same clustered pattern as TPS.

In comparison, our findings show that, like C. camphora and C. kanehirae, C. longepaniculata has similar gene positions and good collinearity on the same chromosomes ( Supplementary Figure S3 ). Tandem repeats, such as MVD1 and MVD2 in C. longepaniculata, as well as MVD3, MVD4, and MVD5 ( Supplementary Figure S3 ), account for the majority of the differences in gene copy numbers between these genomes. Furthermore, structural arrangements lead to gene duplications, as seen in the cases of FPPS2 and FPPS3, which most likely resulted from chromosome segmental duplication events. Notably, an inversion on chromosome 9 reconfigured the TPS gene cluster, resulting in a more densely packed arrangement of these genes, which was also found in the C. camphora genome ( Supplementary Figure S3 ). This data further demonstrates that pathway genes are highly expressed, with at least one gene from each family ( Supplementary Figure S4 ).