Based on the FISH (Fluorescence In Situ Hybridization) results using blue and green light comparison, it was observed that the chromosomes of Michelia champaca and M. alba were fully covered (Figure. S1). This indicates the presence of a close kinship between them. The complete coverage of chromosomes suggests a similarity in the DNA sequence composition and structure between Michelia champaca and M. alba.

The M. alba genome was sequenced using the Illumina HiSeq Xten sequencing platform. K-mer analysis indicated that M. alba has a large genome of approximately 1.8 Gb with 59.1% repetitive elements and is highly heterozygous (4.76%) (Figure. S2). At the same time, the parents were sequenced using the Illumina HiSeq Xten sequencing platform. M. champaca has a large genome of approximately 2.24 Gb but a degree of low heterozygosity (0.38%) (Figure. S3). M. montana has a large genome of approximately 1.47 Gb and has low heterozygosity (0.95%) (Figure. S4). We used the parental Illumina reads to bin the 261 Gb Nanopore reads of the hybrid based on parental origin. A total of 112 G Nanopore reads were assigned to M. champaca, and 130 G Nanopore reads were assigned to M. montana for subsequent haplotype genome assembly (Table S1). A total of 61 G Pacbio Hifi reads were assigned to M. champaca, and 86 G Pacbio Hifi reads were assigned to M. montana for subsequent haplotype genome assembly. The assembled original genome sizes were 2.55 Gb MC (M. champaca) and 2.42 Gb MM (M. montana). respectively. Deduplication was performed by Purge_dups, The obtained genome sizes were 2.23 and 2.19, respectively. Then the raw sequencing data of Pacbio were used for Polish. Using minimap2 + racon strategy to complete three rounds of polish process. The final genome size is 2.19 and 2.13. The contigs and scaffolds of the MC (M. champaca) and MM (M. montana) subgenomes were further scaffolded into 19 chromosomes by Hi-C technology, and the anchored genomes were 2.03 Gb (97.19%) and 2.06 Gb (97.96%), respectively (Figures. S5, S6, Fig. 1a, Table 1). The second-generation genome data of M. alba were mapped to the haplotype genomes MC and MM, with mapping rates of 99.19% and 98.32%, respectively. As a reference, we mapped the second-generation sequencing data of Litchi chinensis, Zea mays, Magnolia biondii, and Liriodendron chinense against the M. alba genome, and the mapping rates were 5.57%, 11.61%, 86.20% and 51.98%, respectively. This indicated a high similarity between M. alba and the MC and MM haplotype genomes. Both assemblies are contiguous (mean Contig N50: 12.425Mb, Scaffold N50:116.01 Mb) and complete (mean BUSCO completeness: 95.5%) (Table S2, Fig. 1c). The corresponding second-generation genome data were compared with the genome data via the Burrows‒Wheeler Aligner (BWA). The mapping rate of the haplotype genome and the Illumina data was 95.68% (MC) and 96.72% (MM). Merqury was used to assess the consensus quality value (QV) and the k-mer completeness of the M. alba genome assembly, which were 34.03 and 96.22%, respectively. These results suggested that the assembled M. alba genome had high completeness and accuracy.

We used the MAKER pipeline and evidence-based and ab initio gene predictions to predict 64,070 protein-coding genes in the MC subgenome assembly and 68,392 protein-coding genes in the MM subgenome assembly. The mean values of gene length, number of exons, and coding sequence length in MC were 1,782.14 bp, 4.6, and 152.929 bp, respectively, and in MM, they were 1,40.95 bp, 4.3, and 148.526 bp, respectively. We also identified 115 and 108 microRNAs (miRNAs), 444 and 492 transfer RNAs (tRNAs), 2,367 and 2,098 small nuclear RNAs (snRNAs), and 174 and 330 ribosomal RNAs (rRNAs) in the subgenomes, respectively. The protein-coding genes in MC and MM had an average gene length of 1,782 bp and 1408 bp, and the average coding DNA sequence (CDS) length was 153 bp and 149 bp, respectively.

The repetitive element percentages in the MC and MM subgenomes were over 74.02% and 71.73%, respectively, comparable to that of Magnolia biondii (Dong et al. 2021) (66.48%) and higher than that of Chimonanthus praecox (Shang et al. 2020a) (~ 45.73%). The percentage of various types of repeats was similar in the haplotype genomes, with the most abundant being Long terminal repeats (LTRs) (34.48% and 30.32%, respectively), followed by long interspersed nuclear elements (LINEs) (0.86% and 1.23%, respectively). The least abundant DNA transposons (0.91% and 0.79%, respectively). (Table S3). Furthermore, using the iTAK tool, we discovered 828 and 859 transcription factors in the MC and MM subgenomes, respectively. (Table S4).

Based on the subgenome, we compared the gene expression differences and structural variations between the two parents and their offspring. The results showed that in the three tissues of petals, stamens and leaves, the number of differential genes in M. champaca and M. alba were 2,702, 1,250 and 2,094, respectively. The number of differentially expressed genes in M. montana and M. alba were 3,322, 2,674, 2,454, respectively. By comparing the second-generation sequencing data of M. alba with the two subgenomes, it can be found that 1,251 insertions and 13,601 deletions were identified between M. champaca and M. alba, while 731 insertions and 13,681 deletions were identified between M. montana and M. alba. It can be seen that the information from the two subgenomes is not significantly different.

The gene families of 6 species, including M. champaca and M. montana, were analyzed. The other four species included were the ancient species Amborella trichopoda (Resource and (TAIR): improved gene annotation and new tools. 2011), Vitis vinifera (Liang et al. 2019), Liriodendron chinense (Chen et al. 2019), and Magnolia biondii. All protein-coding genes from the 6 genomes were clustered into 33,676 gene families (two or more members), of which 5,562 were common to all groups. A total of 1,229 and 881 gene families were found to be specific to M. champaca and M. montana. (Fig. 2a).

We assessed the functional enrichment of specific gene families in the two subgenomes to understand the functional direction of the two subgenomes during speciation. KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis of the MC-specific gene families revealed marked enrichment in genes involved in the biosynthesis of secondary metabolites, starch and sucrose metabolism, glutathione metabolism, biosynthesis of amino acids, fatty acid metabolism and MAPK signaling pathway—plant pathways. On the other hand, MM-specific gene families exhibited a marked enrichment in genes involved in metabolic pathways, carbon metabolism, protein processing in the endoplasmic reticulum, amino sugar and nucleotide sugar metabolism, and glycerophospholipid metabolism (Table S6, Figure. S7).

The position of Magnoliaceae in plant evolution has been extensively studied, but there have been no accurate conclusions. We selected 2 dicot plant species (Vitis vinifera, Arabidopsis thaliana), 3 monocot plant species Elaeis guineensis (Yang et al. 2023), Zea mays (Gui et al. 2020), Oryza sativa (Sakai et al. 2013), 10 Magnoliids species (M. montana, M. champaca, Magnolia biondii, Chimonanthus praecox (Shang et al. 2020b), Cinnamomum kanehirae (Chaw et al. 2019a), Amborella trichopoda (Qin et al. 2021), L. chinense, Piper nigrum (Hu et al. 2019), Magnolia hypoleuca (Zhou et al. 2023) and Magnolia officinalis (Yin et al. 2021), an early angiosperm species (Nymphaea tetragona), an ancient Pteridophyta species (Selaginella tamariscina), and an ancient gymnosperm (Picea abies) for homologous gene analysis to more accurately depict the phylogenetic relationships. After thorough evaluation and selection, we utilized single-copy gene familys to construct the phylogenetic tree. Our findings indicated that monocots and dicots were co-located on one branch, while Magnoliids were in a separate branch of the phylogenetic tree. Furthermore, the phylogenetic tree revealed that M. montana and M. champaca share a close genetic relationship (Table S7, Fig. 2b). In addition, we selected low-copy homologous genes of these 19 species to construct a phylogenetic tree. We reached the same conclusion as the low-copy homologous genes (Figure. S8).

According to our research, the divergence between M. montana and M. champaca occurred approximately 18 million years ago. Additionally, the split between L. chinense and M. montana and M. champaca occurred approximately 122 million years ago. Our investigation also delved into genome evolution, specifically examining whole-genome duplication (WGD) events. We determined the synonymous substitutions per synonymous site (Ks) for various pairs of Magnolia species, including M. champaca vs M. champaca, M. montana vs M. montana, and M. biondii vs M. biondii, as well as comparisons between species such as M. champaca vs M. biondii, M. montana vs M. biondii, and M. champaca vs M. montana. By analyzing the distribution of Ks values, we identified peaks corresponding to genome-wide duplication events. Furthermore, we observed that the Ks distributions of M. champaca and M. montana subgenome paralogues indicated a single WGD event.

To investigate the genome-wide cytosine methylation of the two M. alba haplotypes in various tissues, we analyzed the methylation sites in DNA from both the flowers and leaves of M. champaca and M. montana (as depicted in Fig. 3a). To gain insights into the DNA methylation patterns across different genomic regions, we examined the methylation levels in the gene body (from the transcription start point (TSS) to the transcription termination point (TTS)), as well as the flanking regions 2 kb upstream of the TSS and downstream of the TTS (Fig. 3b and c).

The methylation level of M. champaca flowers and leaves was the highest in the GC sequences, followed by CHH and CHG. Notably, there were few differences in the methylation levels of flowers and leaves. The highest CG methylation levels in flowers and leaves were observed in M. montana, and methylation levels were generally higher in the leaves than in the flowers. The CHH, CG, and CHG methylation sites detected in M. montana flowers were 321,211,379, 317,284,791, and 257,778,276, respectively. The CHH, CG, and CHG methylation sites detected in leaves were 302,713,120, 271,316,538, and 162,325,854, respectively. The CHH, CG, and CHG methylation sites detected in M. champaca flowers were 53,156,789, 52,479,090, and 51,798,809, respectively. The CHH, CG, and CHG methylation sites detected in leaves were 56,094,290, 55,316,685, and 54,198,256, respectively (Table S8). The methylation levels in all three sequence motifs were higher in the leaves of M. montana than in the flowers. In contrast, the opposite was observed in M. champaca, with higher methylation in the flowers.

The regions with significant differences in methylations were studied to further understand M. montana and M. champaca methylation patterns in flowers and leaves. M. montana had 47,381 differentially methylated regions, while M. champaca had 29,907 differentially methylated regions. Moreover, the number of differentially methylated annotated genes was 531 and 343 (Table S9). Thus, M. montana had significantly more differentially methylated genes than M. champaca.

KEGG enrichment analysis of these differentially methylated regions (DMRs) revealed that the M. champaca DMRs were mainly enriched in the biosynthesis of secondary metabolites, biosynthesis of amino acids, and diterpenoid biosynthesis pathways. M. montana DMRs were mainly enriched in the RNA degradation and biosynthesis of amino acid pathways. GO (Gene Ontology) enrichment analysis revealed that M. champaca DMRs were mainly enriched in the carbohydrate metabolic process, response to jasmonic acid, and flavonol biosynthetic process ontologies. M. montana DMRs were mainly enriched in response to wounding, sesquiterpenoid biosynthetic process, and protein phosphorylation ontologies. (Fig. 3e and f, Table S10). Differentially methylated regions were distributed more frequently in the intergenic and promoter regions across the whole genome (Fig. 3d, Table S11).

In this work, we sampled M. champaca and M. montana tissues with yellow and white flower petal colors to elucidate interspecies differences in petal pigments. Using liquid chromatography-electrospray ionization tandem mass spectrometry (LC-ESI-MS/MS), 21 anthocyanin compounds were identified and quantified (against 99 anthocyanin standards) in petals, including 12 delphinidins, 7 pelargonidins, and 2 cyanidins. Eighteen differentially abundant metabolic compounds belonging to anthocyanins were identified between M. champaca and M. montana. Specifically, pelargonidin-3-O-glucoside, pelargonidin 3-sophoroside 5-glucoside, cyanidin 3-xyloside, cyanidin 3-sophoroside-5-glucoside, cyanidin 3-sambubioside, cyanidin 3-glucogalactoside, cyanidin 3-(6' '-succinyl-glucoside), and cyanidin 3-(6' '-malonylglucoside)-5-glucoside were found to be more abundant in M. champaca than in M. montana. One carotenoid compound was also identified but had no significantly different abundance between M. champaca and M. montana. Therefore, anthocyanins represent the most significant differentially abundant metabolites between M. champaca and M. montana. Among these, delphinidins were predominant in M. champaca, with pelargonidin-3-O-glucoside exhibiting the highest concentration (Figure. S9).

Haplotype-based genomic analyses have been frequently employed to understand the relationship between allelic expression and methylation. Based on synteny and annotation, 22,034 pairs were considered reliable allelic genes. Allelic pairs with differential expression in the haplotype genomes were identified using RNA-seq data from two tissues (flower and leaf) of M. champaca and M. montana. A total of 8,099 gene pairs were expressed in the flowers (36.76% of alleles), and 7,095 gene pairs were expressed in the leaves (32.2% of alleles). When differences were identified in the expression between alleles (Table S12), these were considered differentially expressed (Fig. 4a). Statistical analysis of the methylation levels among differentially expressed alleles showed that higher gene expression correlated with a lower methylation level (Sinha et al. 2020) (Fig. 4b), and this relationship was apparent in all tissues (Fig. 4c). It was further verified that the presence of methylated sites inhibits gene expression.

Anthocyanins are pigments that are soluble in water and are frequently and abundantly present in plants. Anthocyanins are responsible for most of the coloring in plant petals. For instance, M. champaca flowers are yellow, M. montana flowers are white, and the flowers of their hybrid offspring, M. alba, are white. In the context of the Magnoliaceae family, it has been discovered that changes in flower color are linked to the abundance of anthocyanins (Liu et al. 2020; Lang et al. 2019). Based on our findings, we compared the methylation levels of anthocyanin synthesis-related genes in M. alba haplotypes. We discovered that M. champaca had higher methylation levels in these genes than M. montana. Additionally, these genes were expressed at lower levels in M. alba flowers and leaves. For example, the homologous genes SSHX090025120.g (MM) and HHX040030052.g (MC) of the key gene 4CL in the flavonoid metabolism pathway, the homologous genes SSHX140010665.g and HHX020009634.g encoding PAL, the homologous genes SSHX100056145.g and HHX040030345.g encoding C4H, and the SSHX090025120.g, SSHX140010665.g, SSHX100056145.g had a much higher expression in flowers and leaves than HHX040030052.g, HHX020009634.g, and HHX040030345.g. However, the methylation levels of these genes in M. montana flowers and leaves were much lower than those in M. champaca (Fig. 4d). From the evidence presented, it is reasonable to speculate that differential methylation plays a role in regulating the expression of flavonoid synthesis genes. This regulation leads to a lower expression level, resulting in the absence of the yellow color observed in the M. montana parental line. This finding serves as a theoretical foundation for further exploration into the mechanisms underlying flower coloration in M. alba.

In addition to their color, the flowers of M. montana and M. champaca differ significantly in appearance. M. montana typically has 9 obovate petals, while M. champaca has 15-20 hard petals with an oblanceolate shape. It is worth noting that the ABCE model is a well-established framework for flower pattern regulation. Four functional genes that are part of the ABCE model were identified in the genome. They were AP1 and AP2 in A, AP3 and PI in B, AG in C, and SEP1, SEP2, SEP3, and SEP4 in E. Based on the gene expression data in petals and stamens in three developmental stages, we observed that the expression of the alleles of these four genes also showed significant differences (Fig. 5a and b, Table S13). In M. montana, the expression of the A genes was very low. AP1 and AP2 were not expressed in the first and third developmental stages of the stamen, and AP2 was not expressed in the petals in the second and third developmental stages. This indicated that A genes played a key role in the process of petal formation. The differential expression of alleles may also be responsible for the flower morphology of M. alba (10 flower petals, lanceolate, petals hard), which is more similar to M. champaca. Furthermore, the analysis of the methylation sites in the promoter regions of the four types of genes revealed that the methylation levels in the CHG sites of M. montana were significantly higher than those of M. champaca. (Fig. 5c).

The M. montana and M. champaca flowers emit a distinctive and delightful aroma. We detected a total of 35 terpenoid compounds in M. champaca and M. montana flowers, including 13 monoterpenes, 14 sesquiterpenes, 3 diterpenes, 1 triterpene, and 4 cyclic ether terpenes. There were significant differences in the composition of these terpenoids between M. champaca and M. montana, as depicted in Figure. S10. This indicates that the distinct aroma profiles of the two species are closely associated with variations in terpenoid composition. In addition, our differential methylation analysis indicated that M. champaca and M. montana had a high abundance of terpenoids, but each had a unique composition of terpenoid compounds. M. champaca had a higher concentration of diterpenes, while M. montana had a higher concentration of sesquiterpenoids. Therefore, we have intensively investigated the key gene family of terpenoid synthesis, terpene synthase (TPS).

We identified TPS gene family members in four species, L. chinense, M. biondii, M. montana, and M. champaca, based on high-quality gene sequence data and annotation. We found that the number of TPS genes in M. montana was significantly higher, with 140 genes present in the genome (Fig. 6a). M. champaca had 121 genes (Fig. 6b), and L. chinense had the lowest number of genes, with only 76 (Table S14, Fig. 6c). We classified the detected TPS genes and discovered the presence of the TPS-c gene in M. champaca. This gene is responsible for diterpene synthesis and is uniquely present in M. champaca compared to M. montana. Additionally, we observed a significant expansion of TPS-a genes in M. montana, which are primarily responsible for synthesizing sesquiterpenoids. By analyzing the transcriptome data from the petals, stamens, and leaves, we observed that gene expression in stamens tended to be higher than that in the other two organs. This suggests that the stamens may play a key role in producing the floral aroma through volatile terpenoid biosynthesis (Fig. 6a and b).

Based on the KEGG and GO enrichment results in M. champaca and M. montana flower and leaf DMR genes, We found that they were enriched in diterpenoid biosynthesis pathways in M. champaca and sesquiterpenoid biosynthesis pathways in M. montana, respectively. (Fig. 6d, Table S15).

Fresh leaves and flowers were collected from the living collections at Hainan University (Haikou, China) and frozen on-site. High-quality genomic DNA was extracted from freshly frozen leaf and flower tissue from one individual M. alba, M. montana, and M. champaca plant using the Plant Genomic DNA Kit (Tiangen) following the manufacturer's instructions.

Fixed plant samples were prepared and treated with specific DNA probes that were labeled with fluorescent dyes or radioisotopes. These labeled probes were then applied to the samples, allowing them to hybridize with complementary target DNA sequences. Following stringent washing steps to remove unbound probes, the samples were examined under a fluorescence microscope or using autoradiography to visualize and detect the hybridization signals.

JELLYFISH 2.1.4 (Marçais and Kingsford 2011) was employed for the K-mer analysis, utilizing K-values of 19, 21, and 25. The analysis focused on the 25-mer oligos, which were fully visualized. The unique sequence, heterozygosity, duplication, and error rates were determined using GenomeScop (Vurture et al. 2017).

A short paired-end Illumina DNA library from the parents M. montana and M. champaca, with a 350 bp insert size (137 × coverage), was sequenced on the Illumina HiSeq 2500 sequencer.

For the M. alba offspring, sequencing with the Nanopore Sequel was performed. Two library-building protocols, the Ultralong reads Protocol for M. alba and a Ligation Library, were used to obtain 75 Gb (13 ×) Nanopore Ultralong reads and 186 Gb (93 ×) Nanopore reads. The official tool Guppy was used for base calling, and a mean_qscore_template value greater than or equal to 7 was used to obtain pass reads. Pass reads can be directly used for the subsequent assembly. The M. alba library produced a total of 261 Gb of data. A PacBio library with an insert size of around 20 kb was then prepared using the SMRTbell Express Template Prep Kit 2.0 from PacBio (Pacific Biosciences, USA). We sequenced the PacBio library on the PacBio Sequel II system (CCS mode), generating 169 Gb clean data (~ 84 ×). A Hi-C library was constructed and sequenced on the Illumina NovaSeq platform for chromosome-level scaffolding, generating 145 Gb of clean data.

We collected samples from the young leaves and stamen tissues of M. alba, M. montana, and M. champaca for transcriptome sequencing. Each tissue sample was isolated and sequenced in three biological replicates to ensure accuracy and reproducibility. To extract the total RNA from the samples, we used the TIANGEN kit with DNASEI, which effectively removes any contaminating DNA. The library was then constructed to prepare it for further processing and analysis. After removing low-quality data, approximately 5.5 Gb of 150-nucleotide paired-end data per sample was used for further RNA-seq analysis.

LC‒MS-based metabolomics analysis was conducted to investigate the metabolic profiles. Six samples of open flowers from M. champaca and M. montana were collected, each in three biological replicates.

For metabolite detection, the sample was first ground into a powder using a grinder. Then, 70% aqueous methanol was added to the powder to extract the metabolites. Metabolites were measured by LC‒MS/MS 6500. Normalization was performed by dividing the relative signal intensity of the metabolites by the intensity of the internal standard (Li et al. 2022).

Canu (Koren et al. 2018) supports using parental short-read sequencing to classify and bin the F1 reads. The long-read progeny sequences were divided into paternal and maternal groups according to haplotype-specific k-mers and assembled separately. Pure third-generation assembly of reads was performed after quality control. The initial genome de novo assembly was done using Hifiasm (Vurture et al. 2017) (v0.15.1) with Nanopore reads, Pacbio Hifi long reads and Hi-C data. To compare the three data generations with the assembled genome, we utilized the minimap2 tool with default parameters. For de-redundancy processing, we employed purge haplotigs (Li et al. 2022) software, utilizing the parameters purge_haplotigs contigcov -l 5; purge_haplotigs purge -a 98. Furthermore, we analyzed the Hi-C libraries using the Juicer pipeline and visualized the results using Juicebox (available at https://github.com/aidenlab/Juicebox).

To predict the gene presence in the existing genome sequence, we performed BUSCO (Koren et al. 2018) prediction using a single copy of the homologous gene in the plant library (eukaryota_odb10).

The prediction of repeat sequences was conducted by aligning the sequences and performing repeat-sequence masking using RepeatMasker (v4.0.6) with the default parameters. The utilized reference repeat libraries involved plant short fragment repeats and annotated repeats of MM and MC. The gene structure prediction was based on a combination of evidence-based and ab initio predictions (denovo). Trinity (v2.4.0) was used to assemble transcripts.

The alignment of RNA-seq data to the reference genome was performed using HISAT2 (v2.0.4) (Cheng et al. 2021), and the assembly was conducted using StringTie (v1.2.2) (Roach et al. 2018). The prediction of gene structures included the utilization of SNAP (Korf 2004), AUGUSTUS (v2.7) (Stanke et al. 2004), and GlimmerHMM (Majoros et al. 2004) as ab initio prediction tools. Maker software and the above results were used to carry out evidence-based assembly.

We utilized protein sequences from MM, MC, and five other angiosperms (A. trichopoda, Liriodendron chinense, Magnolia biondii, and V. vinifera) to identify gene families. This was achieved using OrthoFinder (https://github.com/davidemms/OrthoFinder), employing an all-versus-all BLASTP alignment with an e-value cutoff of 1e-5. Subsequently, gene mapping against the KEGG database allowed us to determine potential gene pathways, while the corresponding Inter-ProScan or Pfam results were used for the extraction of GO terms. Furthermore, OrthoVenn2 was utilized to identify orthologous single-copy genes by using single-copy genes from M. champaca, M. montana, and the aforementioned five angiosperms (A. trichopoda, V. vinifera, L. chinense, and M. biondii).

rRNAs were predicted using RNAmmer (Lagesen et al. 2007) (version 1.2), tRNAs were predicted using tRNAscan-SE (Lowe and Eddy 1997) (version 1.23), and other ncRNA sequences were identified using the Perl program Rfam_scan.pl (version 1.0.4) by inner calling using Infernal (Nawrocki and Eddy 2013) (version 1.1.1). The iTAK (Zheng et al. 2016) program was used to detect known TFs in the MM and MC genomes and the other plants evaluated. The predicted gene set was then used as a query to search the database.

To assess the expansion or contraction of gene families, we utilized CAFE (version 3.0) (Tijl et al. 2006). For further analysis, we selected a minimum of five gene families from various species. Random birth and death patterns were employed to evaluate the changes in gene families across each branch of the phylogenetic tree. The Probabilistic Graphical Model was utilized to determine the probability of gene families transitioning between parent and child nodes during the evolutionary process. This allowed us to analyze the expansion and contraction of gene families across all node types.

The MCScanX (Wang et al. 2012) package was applied to construct syntenic blocks based on well-aligned genes to identify homologous regions between the two haplotypes. We screened the syntenic regions according to the identification rules in the potato diploid genome (Zhou et al. 2020). Many gene pairs may be generated during the identification process, and one gene on the complementary haplotype and its best homolog is generally considered to be an allelic gene pair.

HISAT2 (Kim et al. 2015) was used to construct an index file for the reference genome and mapping. Stringtie (Kovaka et al. 2019), a very commonly used transcriptome expression quantification software, was used to analyze the RNA-seq data and obtain the expression levels in transcripts per million (TPM) of genes on both haplotypes. Two parameters were used to filter the TPM values: (1) Genes showing < 50% difference in expression across three biological replicates were retained. (2) Genes with an average TPM value of 0 for all tissues in three biological replicates were discarded and identified as nonexpressed genes. We compared the differentially expressed allelic genes in two tissues based on genes retained after filtering by TPM value. We considered allelic genes with > twofold expression differences as differentially expressed alleles.

Using tombo (Stoiber et al. 2016) software, the Nanopore sequencing data of the two tissues (flower and leaf) of M. montana and M. champaca were aligned to the reference genome. Then, deepsignal_plant (Ni et al. 2021) software was used to identify methylation sites, and the methylation sites were annotated and plotted. To define the differentially methylated regions (DMRs) between the two haplotypes, we utilized the DSS (Feng et al. 2014) package to obtain the DMRs.

HMMER (Finn et al. 2011) and BLASTP (Altschul et al. 1990) software were used to identify the TPS genes in M. champaca and M. montana. The hmm files were downloaded from the http://pfam.xfam.org/database with serial numbers PF01397 and PF03936. Subsequently, the identified genes were further filtered based on an E value < 0.01. The corresponding Arabidopsis protein data were downloaded from the webpage (https://www.arabidopsis.org/), and the genes were compared based on a 30% identity filter. The intersection of the two was used to obtain the TPS gene database of M. champaca and M. montana. Known TPS sequences were added, muscle (Edgar 2004) software was used to compare the identified sequences, and iqtree (Nguyen et al. 2014) software was used to build tree files and classify TPS genes according to the tree files.

Not applicable.

All authors approve the manuscript and consent to the publication of the work.

The authors declare no competing interests. Xia Zhiqiang is a member of the Editorial Board for Molecular Horticulture. He was not involved in the journal’s review of, and decisions related to, this manuscript.