The wild Dendrobium nobile (voucher specimen: China, Yunnan province, on rock in evergreen broad-leaf forest, alt. 1250 m, 15 April, 2019, GZMU001) plants were collected. The species was identified by comparison with the D. nobile specimen deposited in the herbarium of the National Orchid Conservation Center of China. Green young leaves were frozen for short reads sequencing and PacBio sequencing. Leaf buds from the plants were fixed as described in Wang et al., 2015a and Hi-C library construction is referred to Chen et al. (2020).

For short reads sequencing, genomic DNA was extracted from leaves of the D. nobile plants using a modified cetyltrimethylammonium bromide (CTAB) protocol, and then was shorn by Covaris ultrasonicator. DNA fragments between 300bp and 400bp were then selected by Agencourt AMPure XP-Medium kit for library construction and finally sequencing was performed on the MGISEQ-2000 platform. The raw data was generated from constructed paired-end libraries (PE150), and the clean data was obtained after data filtering, which was carried out by SOAPnuke v1.6.5 software (https://github.com/BGI-flexlab/SOAPnuke) (Chen et al., 2018), with the following parameters: -n 0.02; -l 20; -q 0.4; -Q 2; -i; -G; --seqType 0; –rmdup. The PacBio sequencing was performed on a PacBio Sequel II sequencer by BGI (Shenzhen, China). Furthermore, the data of Hi-C library sequenced by MGISEQ-2000 platform was used for Hi-C analysis.

For assisting gene annotation and gene expression analysis, total RNA was extracted from tender leaves, stems, and roots of three different individuals in the same growing stage of D. nobile using the RNAprep Pure Plant Kit and genomic DNA contamination was removed using RNase-Free DNase I (both from Tiangen), respectively. The integrity of RNA was evaluated on a 1.0% agarose gel stained with ethidium bromide (EB), and its quality and quantity were assessed using a Qubit2.0 Fluorometer and an Agilent 2,100 Bioanalyzer (Agilent Technologies). As the RNA integrity number (RIN) was greater than 7.0 for all samples, they were used in cDNA library construction and Illumina sequencing. The cDNA library was constructed using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB) and 3 μ g RNA per sample, following the manufacturer’s recommendations. The PCR products obtained were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2,100 system. Library preparations were sequenced on the Illumina Novaseq6000 sequencer, generating 150-bp paired-end reads.

Before genome assembly, genome size and heterozygosity were estimated by Jellyfish v.2.1.4 (Marcais and Kingsford, 2011) and GenomeScope (Vurture et al., 2017) based on a 17-K-mer distribution. Then, Canu v 2.2 (Koren et al., 2017) was used to correct the Pacbio raw data and assemble the genome with the following parameters: correctedErrorRate = 0.035 utgOvlErrorRate = 0.065 trimReadsCoverage = 2 trimReadsOverlap = 500 gridOptions = "--mem-per-cpu = 5 g”. Furthermore, pilon v1.22 (--fix bases--mindepth 10 --minqual 20 --diploid) (Walker et al., 2014) was used to correct the assembly using the data generated from the MGISEQ-2000 platform. For chromosome-level assembly, the Hi-C reads were filtered by SOAPnuke software (Chen et al., 2018) with the following parameters: -n 0.02; -l 20; -q 0.4; -Q 2; -i; -G; --seqType 1; –rmdup. Then, the obtained clean reads were compared with the preassembled contigs using Juicer software (Durand et al., 2016). After filtering the results and removing the misaligned reads, 3D-DNA (Dudchenko et al., 2017) software was used to preliminarily cluster, sequence, and direct the pseudochromosomes. Furthermore, Juicer-box was used to adjust, reset, and cluster the pseudochromosomes, and misassemblies and misjoins were manually corrected based on neighboring interactions. For the evaluation of Hi-C assembly results, the final pseudochromosome assemblies were divided into 150 kb bins with equal lengths, and the interaction signals generated by the valid mapped read pairs between each bin were visualized in a heat map. Finally, the completeness and quality of the final assembled genome were evaluated through Benchmarking Universal Single-Copy Ortholog (BUSCO v5.2.2) (Manni et al., 2021) tests.

Repetitive sequences are an important part of the genome and are divided into two types: tandem and interspersed repeats. RepeatMasker v4.0.7 and RepeatProteinMask v4.0.7 software (http://www.repeatmasker.org) were used to identify repetitive sequences based on the RepBase v21.12 database (http://www.girinst.org/repbase). For de novo prediction, a repetitive sequence database was constructed through RepeatModeler (http://www.repeatmasker.org/RepeatModeler/) and LTR_FINDER v1.06 (http://tlife.fudan.edu.cn/ltr_finder/). RepeatMasker software was then used to predict the repeat sequences. Tandem Repeats Finder v4.09 (http://tandem.bu.edu/trf/trf.html) was used to find tandem repeats in the genome. The repeat-masked genome assembly was used for annotating high-quality protein-coding genes through an integration of homology-based, de novo, and transcriptome-based predictions. For homology-based prediction, protein sequences from seven species (Arabidopsis thaliana, Oryza sativa, Asparagus officinalis, A. shenzhenica, Gastrodia elata, P. equestris, and Vanilla planifolia) were used to align the D. nobile genome sequences through Genewise v2.4.1 (Birney et al., 2004) (https://www.ebi.ac.uk/Tools/psa/genewise/). Then, 5,000 complete genes from the homology-based prediction method were used to produce a training model using the Augustus (Stanke et al., 2006) (http://bioinf.uni-greifswald.de/augustus/) and SNAP (Johnson et al., 2008) (http://homepage.mac.com/iankorf/) software. The RNA sequencing data of D. nobile were mapped to the genome sequences through Hisat v2.1.0 and StringTie v1.3.4d (Kim et al., 2015; Pertea et al., 2015). Finally, Maker v2.31.8 (Holt and Yandell, 2011) was used to annotate and integrate the results produced by the three methods. BUSCO v5.2.2 (Manni et al., 2021) was used to evaluate the completeness and quality of the gene models.

Functional annotation of the predicted gene models was carried out by Blast v2.2.26 (Altschul, 1990) software, aligned against the Swissprot (Boeckmann et al., 2003) (http://www.uniprot.org/), TrEMBL (http://www.uniprot.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/), InterPro (Zdobnov and Apweiler, 2001) (https://www.ebi.ac.uk/interpro/), Nr (non-redundant database), KOG (Koonin et al., 2004) (clusters of euKaryotic Orthologous Groups), and Gene Ontology (GO) (Ashburner et al., 2000) databases. For non-coding RNA annotation, tRNAscan-SE 1.3.1 (Lowe and Eddy, 1997) (http://lowelab.ucsc.edu/tRNAscan-SE/) was used to annotate the tRNA sequences. BLASTN was used to search for rRNA. The miRNA and snRNA sequences were predicted by the INFERNAL (Griffiths-Jones, 2004) (http://infernal.janelia.org/) software.

The protein sequences of D. nobile and 17 additional angiosperm species (D. catenatum, D. chrysotoxum, Cymbidium sinense, Amborella trichopoda, Spirodela polyrhiza, A. shenzhenica, A. officinalis, O. sativa, Ananas comosus, Musa acuminata, Phoenix dactylifera, Vitis vinifera, P. equestris, Sorghum bicolor, Populus trichocarpa, A. thaliana, and G. elata) were used for orthologous gene family identification and clustering. BLASTP was used to calculate the similarities between sequence pairs with a cut-off e value of 1e⁻⁵. Finally, OrthoMCL v.2.0.9 (Li et al., 2003) with default parameters was used to cluster the gene families.

The protein sequences of each gene family of 313 identified single-copy gene families were aligned by MUSCLE (Edgar, 2004) (http://www.drive5.com/muscle/). Then, a super-alignment matrix generated from all the alignment results was used for the construction of the phylogenetic tree. Finally, RAxML (Stamatakis, 2014) was used for the construction of the phylogenetic tree of 18 angiosperm species with the GTRGAMMA model; the bootstrap value was 1,000.

PAML MCMCTREE (http://abacus.gene.ucl.ac.uk/software/paml.html) was used to estimate the divergence time. The following constraints were used for time calibrations: i) the B. distachyon and O. sativa divergence time [40.0–54.0 million years (Ma)] (International Brachypodium Initiative, 2010); ii) the divergence time of P. trichocarpa and A. thaliana (100.0–120.0 Ma) (Tuskan et al., 2006); iii) the monocot and eudicot divergence time (lower boundary of 130.0 Ma) (Jaillon et al., 2007); and iv) the time of the earliest-diverging angiosperms (<200.0 Ma) (Magallón et al., 2013).

The Ks distribution was used to infer whole-genome duplication (WGD) events in D. nobile based on paralogous gene pairs and the divergence between species based on orthologues. MCscanX v1.5.2 (Wang et al., 2012) was used to find the collinear regions. BLASTP was used to search for putative paralogous genes in each collinear region within D. nobile and orthologous genes between D. nobile and D. chrysotoxum, D. nobile and V. planifolia, and D. nobile and C. sinense. Furthermore, Codeml, in the PAML package (Yang, 1997) was used with the F3X4 model to calculate the Ks value of each gene pair.

For further analysis of the features of collinear regions, MCscan (https://zenodo.org/record/31631#.XpkUyTOeask) was used to find the collinear regions between D. nobile and D. chrysotoxum, D. nobile and V. planifolia, and D. nobile and C. sinense.

Based on the gene families clustered by OrthoMCL (Li et al., 2003), we filtered the gene families with a number higher than 200 in one species and lower than 2 in other species. CAFÉ software (De Bie et al., 2006) (http://sourceforge.net/projects/cafehahnlab/) was used to determine the expansion and contraction of orthologous gene families combined with divergence times.

All the genes or gene families in the biosynthesis pathways of polysaccharides and alkaloids were identified by HMM or BLASTP searches. The HMM profiles (PF01128.20 for CMS genes, PF01264.22 for CS genes, PF13292.7 for DXS genes, PF02401.19 for HDR genes, PF04551.15 for HDS genes, PF08540.11 and PF01154.18 for HMGS genes, PF00288.27 for PMK genes, PF00275.21 for SHKG genes, PF01202.23 for SK genes, and PF03088.17 for STR genes) were downloaded from Pfam (pfam.xfam.org/), and the query sequences for the BLASTP methods were mostly from homologous genes of Arabidopsis, except for the SKDH genes from Pisum sativum (Weeden, 2018) and DHQS genes from tomatoes (Bischoff et al., 1996). After the candidate homologs were obtained, MAFFT (Katoh and Standley, 2013) and PhyML (Guindon et al., 2005) software were used for sequence alignment and phylogenetic tree construction, respectively.

The HMM profiles for PF01397 (Terpene_synth) and PF03936 (Terpene_synth_C) were downloaded from Pfam (pfam.xfam.org/), and two profiles were used to carry out HMM searches against the protein database for six species (D. nobile, D. chrysotoxum, D. catenatum, P. equestris, A. shenzhenica, and A. thaliana). These sequences were then manually checked, and the sequences with at least one of these domains were retained. The retained amino acid sequences were aligned using MAFFT (Katoh and Standley, 2013). Then, the aligned amino acids were used for phylogenetic tree construction by PhyML (Guindon et al., 2005). The tree was generated by the maximum likelihood method based on the Jones-Taylor-Thornton (JTT) matrix-based model (Jones et al., 1992) and the bootstrap method for phylogeny tests with 1,000 replications. The HMM profiles for PF00067.23 (CYP450) were downloaded from Pfam (pfam.xfam.org/). HMM searches were used to obtain the homologs in the protein database for seven species (D. nobile, D. chrysotoxum, D. catenatum, P. equestris, A. shenzhenica, A. thaliana, and O. sativa). The next steps were the same as those for TPS gene family identification, except for the rice genome. The identified CYP450 genes in A. thaliana and O. sativa in this study were further confirmed with the Arabidopsis Cytochrome P450 database (http://www.p450.kvl.dk/At_cyps/family.shtml) and rice (Wei and Chen, 2018). All gene expression analyses were carried out using Salmon v1.3.0 (Patro, et al., 2017) with the default settings.

D. nobile (Figure 1) has a karyotype of 2N = 2X = 38 (Zheng et al., 2018). For estimating the D. nobile genome size, a total of 130.62 Gb of raw data with 300-400 bp insert libraries were generated by MGISEQ-2000 sequencing and 122.52 Gb of clean data was obtained after data filtering (Supplementary Table S1). The estimated genome size was 1.16 Gb, with 1.35% heterozygosity based on K-mer analysis (Supplementary Figure S1). To obtain the assembly, 96.91 Gb (coverage of 83.54×) of PacBio sequencing data and 5.19 million subreads with an N50 read length of 21.4 kb was generated (Supplementary Table S1). Based on the PacBio and MGISEQ-2000 sequencing data, we then used 245.48 Gb of raw data from a Hi-C library (Supplementary Table S1) to reconstruct physical maps by recoding and clustering the assembled scaffolds into 19 pseudochromosomes, which represented the 19 chromosomes in the haploid genome of D. nobile (Figure 2; Supplementary Table S2). Finally, 211.34 Gb of clean data was produced using the data filtering. After the Hi-C analysis, 1,192,377,146 bp sequences were mapped to 19 pseudochromosomes, accounting for 99.45% of the raw assembly (1.19 Gb). The length of the 19 pseudochromosomes ranged from 37.78 to 95.36 Mb, with a scaffold N50 value of 64.46 Mb (Table 1; Supplementary Tables S2, S3). Furthermore, BUSCO estimation indicated that the completeness of the gene set of the assembled genome was 96.22%, while mapping short reads back to the assembly indicated that the completeness was 98.82% (Supplementary Tables S4, S5). The chromatin interaction data showed a high level of chromatin interaction between linked sequences and low chromatin interaction between non-linked sequences, suggesting the high quality of the Hi-C assembly (Figure 3).

A total of 29,476 protein-coding genes were annotated in D. nobile (Supplementary Table S6). Among them, 98.28% of ≥50% CD regions overlapped based on homolog or de novo data (Supplementary Table S7), with 93.43% exhibiting high completeness (Supplementary Table S8).

The average lengths of genes and introns in D. nobile were similar to those in other orchid species, except for A. shenzhenica (Zhang et al., 2017), and much higher than those in most other angiosperms (Zhang et al., 2017; Supplementary Figure S2). This suggests that the gene length is positively correlated with the intron length, and that genes and introns in orchids except for A. shenzhenica have a longer average length, which may be a unique characteristic. Furthermore, 76 microRNAs, 386 transfer RNAs, 958 ribosomal RNAs, and 457 small nuclear RNAs were identified in the D. nobile genome (Supplementary Table S9).

In this study, it was estimated that 61.07% of the D. nobile genome consisted of repetitive sequences (Supplementary Table S10), similar to the 62% in P. equestris and 62.81% in D. chrysotoxum, and lower than the 78.1% in D. catenatum (Cai et al., 2015; Zhang et al., 2016). Furthermore, the notably higher percentage of repeats predicted by the de novo method compared to the homology-based search method indicates that D. nobile has many characteristic repeats (Supplementary Table S10). Among these elements, long terminal repeats (LTRs) accounted for 51.31% of the genome (Supplementary Table S11), similar to the 53.15% in D. chrysotoxum and higher than the 46% in D. catenatum, 46.47% in P. equestris, and 22.06% in A. shenzhenica (Zhang et al., 2016; Zhang et al., 2017).

In addition, 27,765 (94.20%) predicted genes were functionally annotated (Supplementary Table S12). Among them, 27.601 (93.64%) and 25,870 (87.77%) genes were annotated to the Nr and TrEMBL databases, respectively, (Supplementary Table S14). The number of annotated genes was 24,044 (81.57%), 20,215 (68.58%), and 19,855 (67.36%) in the Interpro, KEGG, and Swissprot databases, respectively (Supplementary Table S12).

For the phylogenetic relationship and divergence times among different plant species, a high-confidence phylogenetic tree and the estimated divergence times of 18 different plant species based on genes extracted from a total of 313 single-copy families were constructed (Supplementary Figures S3, S4; Supplementary Table S13). As expected, D. nobile, D. catenatum, and D. chrysotoxum were sisters to P. equestris, forming an Epidenroideae clade. G. elata and A. shenzhenica were located at the base of the Orchidaceae branches (Supplementary Figure S4). The estimated Orchidaceae divergence time was 113.3 (104.3–121.0) Mya, the divergence time of the subfamily Apostasioideae was 76.4 (61.1–90.8) Mya, and the divergence time between D. nobile and P. equestris was 41.7 (29.2–54.9) Mya (Figure 4). Finally, D. chrysotoxum appeared at 11.9 (6.5–18.1) Mya, and the divergence time between D. nobile and D. catenatum was 5.6 (2.9–8.6) Mya. Then, the expansion and contraction of orthologous gene families were determined. According to the results (Figure 5), 110 gene families were expanded in the lineage leading to the Orchidaceae, whereas 1,566 gene families were contracted. Furthermore, 595 gene families expanded and 321 gene families contracted, leading to the Dendrobium genus. In D. nobile, 859 gene families were expanded compared with 417 in D. catenatum, 721 in D. chrysotoxum, and 872 in P. equestris. For the contracted gene families, 325 gene families were contracted in D. nobile compared with 1,010 in D. catenatum, 1,510 in D. chrysotoxum, and 900 in P. equestris. In the D. nobile clade, 859 gene families were expanded, including 3,134 genes, and 325 gene families were contracted, including 213 genes. To further investigate the functions of the expanded gene families and their species specificity, KEGG enrichment analysis was conducted for the expanded gene families. The KEGG term “Biosynthesis of other secondary metabolites” was found to be the most significantly enriched, and contained the most genes (Supplementary Table S14). This may be related to the synthesis of specific medicinal ingredients in D. nobile.

To detect the occurrence of WGDs in D. nobile, the Ks distribution pattern of paralogous genes in collinear regions was analyzed. There were two peaks, and the Ks values were 0.8 and 1.5 (Figure 6), suggesting the occurrence of two recent WGD events in D. nobile. To determine the time of the two WGD events, the Ks distribution of the homologous genes between D. nobile and P. equestris, A. shenzhenica, D. chrysotoxum, and G. elata was further analyzed. The most recent WGD event in D. nobile (Ks = 0.8) occurred before the divergence between D. nobile and A. shenzhenica, suggesting that the WGD events were shared among all extant orchids (Zhang et al., 2017). The second recent WGD event in D. nobile (Ks = 1.5) was shared with most monocots (Zhang et al., 2017). For the collinear regions, the syntenic relationships among D. nobile, V. planifolia, C. sinense, and D. chrysotoxum were further analyzed. There were 25,166 collinear gene pairs located in 1,404 collinear regions between D. nobile and V. planifolia. 31,498 collinear gene pairs located in 120 collinear regions between D. nobile and D. chrysotoxum, and 30,650 collinear gene pairs located in 430 collinear regions between D. nobile and C. sinense (Supplementary Figures S5–S7; Supplementary Table S15).

The key genes encoding enzymes in the potential polysaccharide synthesis pathway were identified in D. nobile, D. catenatum, D. chrysotoxum, A. shenzhenica, P. equestris, and A. thaliana (Figure 7; Supplementary Table S16). All six gene families were multigene families. There were seven alkaline/neutral invertase (NI) genes, two PGM genes, three SPS genes, four SUS genes, two UGP genes, and five UGE genes in D. nobile. Based on the phylogenetic analysis, these gene families could be further divided into several branches. The number of genes in different branches among these species varied greatly, especially those in the NI, PGM, and UGE gene families (Figure 7, in red font). The number of NI genes in the cytosolic branch, mitochondria branch, and chloroplastic branch of D. nobile was five, two, and one, respectively, compared with seven, one, and one in D. chrysotoxum, respectively, and three, three, and one in D. catenatum, respectively (Supplementary Table S17). For PGM genes, one gene in D. nobile was clustered into the pPGM branch, and one gene were divided into the cPGM branch (Supplementary Table S16). Interestingly, there were five PGM genes of D. chrysotoxum on the cPGM branch and one gene on the pPGM branch (Supplementary Table S16). There were five and zero members divided into group 1 and group 2 in the UGE gene family, respectively, while there were one to two different values among other species (Supplementary Table S16). The gene number variations in different branches of each gene family may be related to the different amounts and constituent polysaccharides among Dendrobium species. Furthermore, for the SPS4F branch in the SPS gene family, three homologous genes were identified in three Dendrobium species, and no genes were found in other orchid species (A. shenzhenica and P. equestris), suggesting a Dendrobium-specific gene (Figure 7, in green font; Supplementary Table S16). The gene expression patterns in different tissues of the three Dendrobium species were also analyzed, but most were not significantly different among tissues.

To determine the genetic variation of the alkaloid biosynthesis pathway, all key genes (20 genes or gene families) upstream of the alkaloid biosynthesis pathway were identified among three Dendrobium species (D. nobile, D. catenatum, and D. chrysotoxum) and three other species (A. shenzhenica, P. equestris, and A. thaliana) (Supplementary Table S17)). Most were single-copy genes, except for the DHS, DXS, and HMGR gene families. Interestingly, no SKDH, HMGR, or MVD genes were found in D. chrysotoxum, which hinted at a significant variation in genes or the possibility of alternative pathways playing a dominant role in this process. The gene expression patterns in different tissues among the three Dendrobium species were also determined. Three genes or gene families in the shikimate pathway exhibited higher expression in stems than in leaves: Dnobile10G00647.1 (DHS), Dnobile17G01430.1 (DHS), Dnobile08G00887.1 (DHQS), and Dnobile02G00904.1 (DHD-SKDH) (Figure 8, red five-pointed star; Supplementary Figure S8). The Dnobile08G00649.1 (HMGS) and Dnobile01G00682.1 (MVD) genes in the MVA pathway also had higher expression in stems (Figure 8, red five-pointed star). All genes in the MEP pathway were expressed equally in stems and leaves. For D. chrysotoxum, most genes had higher expression in leaves than in stems (Supplementary Figure S9), such as Guchui-Maker62728 (DHS), Guchui-Maker69398 (DXR), Guchui-Maker111118 (MCS), Guchui-Maker68853 (HDS), and Guchui-Maker109159 (HDR). For STR genes, Dnobile14G00298.1 (STR12/13) was mainly expressed in stems (Supplementary Figure S10), while other members mainly low or equally expression between the two tissues. Interestingly, Guchui-Maker98640 (STR9) was mainly expressed in stems, and Guchui-Maker96976 (STR10) was mainly expressed in leaves (Supplementary Figure S10). Dendrobine is mainly produced in the stems of D. nobile (Li et al., 2019) and is rare in D. chrysotoxum. These opposite gene expression patterns may be related to the difference between D. nobile and D. chrysotoxum in terms of alkaloid (dendrobine) biosynthesis.

Dendrobine, a sesquiterpenoid alkaloid, is the dominant type of alkaloid in D. nobile. TPSs catalyze the biosynthesis of sesquiterpenes (C15) using farnesyl diphosphate (FPP) as a substrate (McGarvey and Croteau, 1995). TPSs may have originated from isoprenyl diphosphate synthase genes, which are involved in dendrobine biosynthesis (Jiang et al., 2019; Wang et al., 2020). CrGES (a gene encoding geraniol synthase that is present in Catharanthus roseus, and is the homologue of TPS02) and TPS21 genes are involved in dendrobine biosynthesis (Kuang et al., 2016; Wang et al., 2020).

In this study, 51 TPS genes were identified in D. nobile, which was more than the 48 identified in D. chrysotoxum and the 42 identified in D. catenatum (Figure 9). The number of TPS genes in Dendrobium was much higher than that in A. shenzhenica and P. equestris, which indicated a relationship with dendrobine biosynthesis. The TPS gene family was divided into five subfamilies, TPS-a, TPS-b, TPS-c, TPS-e/f, and TPS-g, and the TPS-b subfamily was further clustered into TPS-b-I, TPS-b-II, and TPS-b-III. The TPS-a and TPS-b subfamilies showed different evolutionary patterns between monocots and dicots. Interestingly, there were 21, 16, and 17 TPS-a genes in D. nobile, D. chrysotoxum, and D. catenatum, respectively. As the TPS-a subfamily encodes only sesqui-TPSs (Jiang et al., 2019), a higher number of TPS-a genes in D. nobile may be an important factor in the higher production of dendrobine compared with other Dendrobium species (Li et al., 2019). Based on gene expression pattern analysis, genes of the TPS-a subfamily were mainly expressed in the stem (Supplementary Figure S11), which further confirmed that these genes contributed to the production of more dendrobine in D. nobile through the existence of more genes with higher expression in the stems. Furthermore, the gene numbers of TPS-b-II and TPS-b-III were three and 13 in D. nobile, respectively, compared with four and 10 in D. catenatum, respectively, while the numbers were 14 and seven in D. chrysotoxum, respectively. There were several TPS genes with higher expression in the leaves than in the stem (Supplementary Figure S11), such as Dnobile18G01782.1, Dnobile18G01781.1, Dnobile18G01775.1, and Dnobile18G01780.1. As the TPS-b subfamily encoding monoTPSs (Jiang et al., 2019), the opposite distribution and gene expression pattern may have contributed to the differences in the amounts and constitution of alkaloids between D. nobile and D. chrysotoxum in different tissues.

Cytochrome P450 monooxygenases (CYP450s) are integral components in terpenoid and alkaloid pathways (Schuler, 2015). Biosynthetic CYP450s in these pathways can be considered organism-specific (Schuler, 2015). In this study, all the CYP450 genes were identified in three Dendrobium species (D. nobile, D. chrysotoxum, and D. catenatum), two other orchid species (P. equestris and A. shenzhenica), A. thaliana, and rice. There were 123, 170, 228, 257, and 210 CYP450 gene members in A. shenzhenica, P. equestris, D. nobile, D. chrysotoxum, and D. catenatum, respectively. The number of genes in Dendrobium species was much higher than that in P. equestris and A. shenzhenica, mainly because of tandem gene duplication (Supplementary Figure S12).

CYP450 families are grouped into the 10 plant CYP450 clans (Nelson et al., 2004), and many duplications and divergence events, such as chemical defense pathways in particular plants or groups of plants, occur mainly in the expanded multiple-family clans. The species-specific synthesis and catabolism pathways of terpenoids are mainly caused by gene neofunctionalizations within the CYP85 clan (Schuler, 2015). In this study, the differences in CYP450 genes among species were also mainly detected in the CYP85 clan. The CYP720 subfamily, which mediates the oxygenation of monoterpenes (myrcene and pinenes), sesquiterpenes (farnesene), and diterpenes (abietadienol and abietic acid) (Ro et al., 2005; Zulak and Bohlmann, 2010; Hamberger et al., 2011), was only identified in dicots in a previous study (Wei and Chen, 2018). In this study, the orthologous genes of CYP720 in orchids were identified (Supplementary Figure S12, red). Differences among CYP85 clans may involve the dendrobine signaling pathway in D. nobile.

Furthermore, other differences in orchid-specific related traits in CYP450 genes were also identified. CYP78 is involved in regulating organ size and cell proliferation (Nagasawa et al., 2013; Wang X. et al., 2015). More genes expanded by tandem duplication were found in Dendrobium species (Supplementary Figure S12, pink), such as Dnobile03G00500.1, Dnobile03G00501.1, and Dnobile03G00502.1 in D. nobile, GuChui-Maker69281, GuChui-Maker69269, GuChui-Maker69278, GuChui-Maker23937, and GuChui-Maker23933 in D. chrysotoxum, which may be involved in the formation of large capsules and swollen stems. CYP715, as a single gene family, plays a role as a key regulator in flower maturation, synchronizing petal expansion, and volatile emission in Arabidopsis (Liu et al., 2015). Interestingly, there were two to three CYP715 homologous genes in three Dendrobium species and P. equestris (Supplementary Figure S12, blue), indicating that gene expansion occurred in the Epidendroideae subfamily, potentially as a result of floral development.