Tripterygium wilfordii plants were collected from the experimental fields of Shanghai Chenshan Botanical Garden (31° 04′ 30.00′′ N, 121° 10′ 58.93′′ E) and cultured in a greenhouse by cutting propagation. All materials used for genome sequencing originated from a single plant (grown in the greenhouse of our laboratory, voucher TW1). For RNA sequencing (RNA-seq), tissues from roots (R), stems (S), young leaves (YL), mature leaves (L), flower buds (FB) and flowers (F) were harvested with three independent biological replicates.

Total DNA was isolated from leaves using the modified cetyltrimethylammonium bromide (CTAB) method [23]. DNA purity was checked by electrophoretic analysis on a 1% agarose gel and using a NanoPhotometer spectrophotometer (IMPLEN, CA, USA). The DNA concentration was determined using a Qubit 2.0 fluorometer (Life Technologies, CA, USA).

For Illumina sequencing, qualified DNA was fragmented using a Covaris device (MA, USA). Fragmented DNA was end-repaired; poly(A) tail and adaptor addition was performed using the Next Ultra DNA Library Prep Kit (NEB, MA, USA), then the appropriate samples were selected by electrophoretic analysis. The size-selected product was PCR-amplified, and the final product was purified and validated using AMPure XP beads (Beckman Coulter, CA, USA) and an Agilent Bioanalyzer 2100. Using the HiSeq 2500 platform, 150-bp paired-were sequenced. Clean data were obtained by removing adaptor reads, unidentified nucleotides (N) and low-quality reads from the raw reads, and Q20, Q30 and GC content of the clean data were calculated for quality assessment (Table 1). We estimated the genome size by performing k-mer frequency analysis. The k-mer frequencies (k-mer size = 17) were obtained using Jellyfish v2.2.7 [24] with jellyfish count -G 2 -s 5G -m 17 and jellyfish stats as the default parameters.

For long-read sequencing, qualified DNA was sheared into fragments in a g-TUBE (Covaris, MA, USA) by centrifugation, and quantity and quality were controlled by an Agilent Bioanalyzer 2100. To construct a sequencing library, the fragmented DNA was end-repaired and poly(A) tail and adaptor addition was performed using Next Ultra II End Repair/dA-Tailing Module, Next FFPE DNA Repair Mix and Next Quick Ligation Module (NEB, MA, USA), respectively, according to the manufacturer’s instructions. The final product was validated using an Agilent Bioanalyzer 2100. Finally, the qualified DNA library was sequenced using Oxford Nanopore Technology (ONT) on the PromethION platform.

De novo genome assembly was carried out using NextDenovo v2.3.0 [25]. The correct_option parameters used were: read_cutoff = 1k, seed_cutoff = 28087, pa_correction = 20, seed_cutfiles = 100, sort_options = -m 15g -t 8 -k 40, minimap2_options_raw = -x ava-ont -t 8. The assemble_option parameters used were: random_round = 20, minimap2_options_cns = -x ava-ont -t 8 -k17 -w17, nextgraph_options = -a 1.

Racon v1.3.1 [26] and Pilon v1.22 (Pilon, RRID:SCR_014731) [27] were used for error correction with ONT data and Illumina data, respectively. Error correction was performed three times with default parameters. The completeness of the genome assembly was assessed using Benchmarking Universal Single-Copy Orthologs (BUSCO) v3.0.2 (BUSCO, RRID:SCR_015008) [28] with the parameters: -m genome -c 15 -sp arabidopsis. Assembly accuracy was evaluated using Burrows–Wheeler Aligner (BWA) software (version: 0.7.8-r455) [29] to align Illumina reads back to the genome. Variant calling was performed using SAMtools (version: 0.1.19-44428cd, SAMTOOLS, RRID:SCR_002105) [30, 31] with parameters: -m 2 -F 0.002 -d 110 -u -f. To assess the genome assembly quality, transcriptome data were assembled using Trinity (Trinity, RRID:SCR_013048) [32], then mapped back to the scaffolds using BLAT (BLAT, RRID:SCR_011919) [33].

For Bionano sequencing, genomic DNA (molecules > 300 kb) of leaves from living T. wilfordii plants was extracted using the Plant DNA Isolation Kit (Bionano Genomics, CA, USA). Using the NLRS DNA Labeling Kit (Bionano Genomics, CA, USA), DNA molecules were digested with Nt.BspQI endonucleases (determined after evaluation by electronic digestion), and fluorescently labeled. Labeled DNA molecules were electrophoretically stretched into linearization by Saphyr Chip (Bionano Genomics, CA, USA), passed through the NanoChannels [34], and then captured on the Saphyr platform with a high-resolution camera. Raw image data were first converted to digital representations of the motif-specific label pattern, then analyzed using Bionano Solve v3.1 [35] and its in-house scripts. Bionano data were compared with the draft genome (Nanopore version) with the parameters: -U, -d, -T, 3, -j, 3, -N, 20, -I, 3, and scaffolds were generated by connecting contigs with the parameters: -f, -B, 1, -N, 1.

Hi–C library preparation and sequencing were based on a protocol described previously, with some modifications [36–38]. Leaves from living T. wilfordii plants were treated with 1% formaldehyde solution to fix chromatin. Approximately 2 g of fixed tissue was homogenized with liquid nitrogen, resuspended in nucleus isolation buffer and filtered with a 40-nm cell strainer. Extracted chromatin was cut with the HindIII restriction enzyme (NEB, MA, USA), end-filled, then labeled with biotin. After ligation with T4 DNA ligase (NEB, CA, USA) and reversal of crosslinking by proteinase K, DNA was purified, cleaved into 350-bp fragments and end-repaired. DNA fragments labeled with biotin were separated using Dynabeads M-280 Streptavidin (Thermo Fisher, MA, USA), purified, and end-repaired. A-tails were added and adaptors were ligated, and the sequences were amplified by PCR to generate Hi–C libraries. Finally, the qualified libraries were sequenced on an Illumina platform. Clean data were obtained by removing adaptor reads, unidentified nucleotides (N) and low-quality reads from the raw reads, and Q20, Q30 and GC content of the clean data were calculated for quality assessment (Table 2). Clean data were first mapped to the draft genome using BWA software (version: 0.7.8-r455) [29]. After removal of PCR duplicates and unmapped reads using SAMtools v1.9 [39], based on the numbers of interacting read pairs, contigs were clustered and ordered into chromosome groups using LACHESIS (version 201701) [40] with the parameters: RE_SITE_SEQ = GATC, CLUSTER_N = 23, CLUSTER_CONTIGS_WITH_CENS =  −1, CLUSTER_MIN_RE_SITES = 388, CLUSTER_MAX_LINK_DENSITY = 3, CLUSTER_NONINFORMATIVE_RATIO = 0, CLUSTER_DRAW_HEATMAP = 1, and CLUSTER_DRAW_DOTPLOT = 1.

Total RNA was extracted from collected tissues using the RNAprep Pure Plant Kit (TIANGEN, Beijing, China). Qualified RNA from each sample was used to generate sequencing libraries using the NEBNext Ultra RNA Library Prep Kit for Illumina sequencing (NEB, MA, USA), following the manufacturer’s instructions. mRNA was purified from total RNA using poly(T) oligo-attached magnetic beads, then cleaved into short fragments that were used as templates for cDNA synthesis. After purification, repair, adenylation and adaptor ligation of the 3′ end, 150 to 200-bp cDNA fragments were separated for PCR amplification. Finally, libraries were quality controlled using an Agilent 2100 Bioanalyzer and qPCR, then sequenced on the Illumina HiSeq 2500 platform. Raw reads with adapter, poly(N) and low-quality reads were removed to generate clean data, and Q20, Q30 and GC content of the clean data were calculated for quality assessment (Table 3). RNA-seq data was mapped back to the genome assembly of T. wilfordii using HISAT v2.0.4 with default parameters [41]. The read numbers of each gene were counted using HTSeq v0.6.1 with the parameter: –m union [42]. Fragments per kilobase of transcript per million fragments mapped (FPKM) of each gene was calculated based on the length of the gene and number of read counts mapped to this gene [43].

For full-length transcriptome sequencing by PacBio, the best quality RNA samples of each tissue were mixed together to build an isoform sequencing library using the Clontech SMARTer PCR cDNA Synthesis Kit and the BluePippin Size Selection System protocol, as described by Pacific Biosciences (PN 100-092-800-03). Samples were then sequenced on the PacBio Sequel platform. Sequence data were processed using SMRTlink 7.0 software [44] with the parameters: –minLength 50, –maxLength 15000, –minPasses 1. Error correction was achieved using the Illumina RNA-seq data with LoRDEC v0.7, with the parameters: -k 23, -s 3 [45]. Redundancy in the corrected consensus reads was removed by CD-HIT v4.6.8 [46], with the parameters: -c 0.95, -T 6, -G 0, -aL 0.00, -aS 0.99, -AS 30 to obtain the final transcripts for the subsequent analysis.

Homolog alignment and de novo prediction were applied for repeat annotation. For homolog alignment, the Repbase database employing RepeatMasker software v4.0.7 (RepeatMasker, RRID:SCR_012954) and its in-house scripts (RepeatProteinMask v4.0.7) was used with default parameters to extract repeat sequences [47]. For de novo prediction, LTR_FINDER v1.0.7 (LTR_Finder, RRID:SCR_015247) [48], RepeatScout v1.0.5 (RepeatScout, RRID:SCR_014653) [49], and RepeatModeler v1.0.3 (RepeatModeler, RRID:SCR_015027) [50] were used with default parameters to build a de novo repetitive element database for repeat identification. Tandem repeats were also extracted by de novo prediction using TRF v4.0.9 [51].

A combined strategy based on homology, gene prediction, RNA-seq and PacBio data was used to annotate gene structure. For homolog prediction, sequences of proteins from six species, including Arabidopsis thaliana, Vitis vinifera, Medicago truncatula, Cucumis sativus, Ricinus communis, and Glycyrrhiza uralensis, were downloaded from Ensembl/National Center for Biotechnology Information (NCBI)/DNA Database of Japan (DDBJ). Protein sequences were aligned to the genome using TblastN v2.2.26 [52] (E-value ≤ 1 × 10−5), then the matching proteins were aligned to homologous genome sequences for accurate spliced alignments with GeneWise v2.4.1 (GeneWise, RRID:SCR_015054) [53]. De novo gene structure identification was based on Augustus v3.2.3 (Augustus, RRID:SCR_008417) [54], GlimmerHMM v3.04 (GlimmerHMM, RRID:SCR_002654) [55], and SNAP (2013-11-29) [56]. Based on the above prediction results, RNA-seq reads from different tissues, and PacBio reads, were aligned to the genome using HISAT v2.0.4 (HiSat2, RRID:SCR_015530) [41] and TopHat v2.0.12 (TopHat, RRID:SCR_013035) [57] with default parameters to identify exon regions and splice positions. Alignment results were then used as input for Stringtie v1.3.3 (StringTie, RRID:SCR_016323) [58] with default parameters for genome-based transcript assembly. Alignment results were then integrated into a nonredundant gene set using EVidenceModeler v1.1.1 and further corrected with Program to Assemble Spliced Alignment (PASA) to predict untranslated regions and alternative splicing to generate the final gene set [59].

According to the final gene set, gene function was predicted by aligning the protein sequences to Swiss-Prot [60] and the Non-Redundant Protein Sequence Database (NR) (version 20190709) [61]. The motifs and domains were annotated using InterProScan70 v5.31 by searching against the Protein Families Database (Pfam) [62], Kyoto Encyclopedia of Genes and Genomes (KEGG, version 20190601) [63], and Integrative Protein Signature Database (InterPro) v32.0 [64] using Blastp (E-value ≤ 1 × 10−5). Gene Ontology (GO) IDs for each gene were assigned according to the corresponding InterPro entry.

Noncoding RNA was annotated using tRNAscan-SE v1.4 (for tRNA) [65] or INFERNAL v1.1.2 with default parameters (for miRNA and snRNA) [66]. rRNA was predicted by BLAST using rRNA sequences from A. thaliana and O. sativa as references, which are highly conserved among plants.

Gene family clustering of 12 species, including T. wilfordii, A. thaliana, Citrus sinensis, V. vinifera, Glycine max, M. truncatula, G. uralensis, C. sativus, Populus trichocarpa, R. communis, Oryza sativa, and Amborella trichopoda, was inferred through all-against-all protein sequence similarity searches using OthoMCL v1.4 [67], with the parameters: -mode 3 and -inflation 1.5. Proteins containing fewer than 50 amino acids were removed, and only the longest predicted transcript per locus was retained.

Single-copy orthologous genes were retrieved from the 12 species and aligned using MUSCLE v3.8.31 with default parameters [68]. All alignments were combined to produce a super-alignment matrix, which was used to construct a maximum likelihood (ML) phylogenetic tree using RAxML v8.2.12 [69] with the parameters: cds: -m GTRGAMMA -p 12345 -x 12345 -#100 -f ad -T 20, pep: -m PROTGAMMAAUTO -p 12345 -x 12345 -#100 -f ad -T 20.

Divergence times between species were calculated using the MCMCtree v4.9 program implemented for phylogenetic analysis by maximum likelihood (PAML) with the default parameters [70]. The following calibration points were applied: M. truncatula–G. uralensis (15–91 million years ago, Mya), G. max–M. truncatula (46–109 Mya), G. max–C. sativus (95–135 Mya), A. thaliana–C. sinensis (96–104 Mya), P. trichocarpa–R. communis (70–86 Mya), A. thaliana–P. trichocarpa (98–117 Mya), C. sativus–R. communis (101–131 Mya), V. vinifera–A. thaliana (107–135 Mya), V. vinifera–O. sativa (115–308 Mya), and O. sativa–A. trichopoda (173–199 Mya). These calibrations were extracted from TimeTree [71].

Expansion and contraction of gene families were analyzed by using CAFÉ v4.2 [72] with the parameters: -p 0.05 -t 4 -r 10000. To avoid false positives, results were filtered and the enrichment results screened with a family-wide P-value < 0.05 and Viterbi P-values < 0.05.

The hidden Markov model (HMM) profile of Pfam PF06200 [62] was used to extract full-length CYP candidates from the T. wilfordii genome by the HMM algorithm (HMMER) [71], filtering by a length between 400 and 600 amino acids [74].

Multiple sequence alignments and phylogenetic tree construction were performed using MEGA X [75], with either the neighbor-joining or ML method with a bootstrap test (n = 1000 replications).

Gene expression pattern analysis was performed using Short Time-series Expression Miner (STEM) software [76] on the OmicShare tools platform [77]. The parameters were set as follows: the maximum unit change in model profiles between time points was 1; the maximum output profile number was 20 (similar profiles were merged); the minimum ratio of fold change of differentially expressed genes (DEGs) was no less than 2.0, and the P-value was <0.05.

The complete open reading frames (ORFs) of the putative CYP genes were amplified using the primers listed in Table 4, with cDNA from T. wilfordii root used as the template. According to the manufacturer’s instructions, fragments were cloned into the entry vector pDONR207 and yeast expression vector pYesdest52 using the Gateway BP Clonase II Enzyme Kit and LR Clonase II Enzyme Kit (Invitrogen, MA, USA), respectively.

Friedelin, 29-hydroxy-friedelan-3-one, and celastrol were purchased from Yuanye-Biotech (Shanghai, China), and polpunonic acid and wilforic acid A were purchased from Weikeqi-Biotech (Sichuan, China). Friedelin was dissolved in dimethyl sulfoxide (DMSO)/isopropanol (v/v = 1:2) following 30 min of ultrasonication in a water bath, while 29-hydroxy-friedelan-3-one, celastrol, polpunonic acid and wilforic acid A were dissolved in methanol.

Plant tissue was ground into powder in liquid nitrogen then freeze dried. Fifty milligrams of sample was suspended in 2 mL of 80% (v/v) methanol, set overnight at room temperature, then extracted in an ultrasonic water bath for 60 min. After centrifugation at 12,000g for 2 min, the supernatant was filtered through a 0.2-μm Millipore filter before liquid chromatography–mass spectrometry (LC–MS) analysis.

Levels of celastrol and wilforic acid A were analyzed using an Agilent 1260LC-6400 QQQ (triple quadrupole mass spectrometer). Chromatographic separation was carried out on an Agilent Eclipse XDB-C18 analytical column (4.6 × 250 mm, 5 μm) with a guard column. The flow rate of the mobile phase consisting of 0.1% (v/v) formic acid in water (A) and acetonitrile (B) was set to 0.8 mL/min. The gradient program was as follows: 0–12 min, 10–60% B; 12–17 min, 70% B; 17–25 min, 95% B; 25–28 min, 95% B; 28–29 min, 5% B; 29–35 min, 5% B. The detection wavelength of celastrol was 425 nm, and UV spectra from 190–500 nm were also recorded. The injection volume was 10 μl and the column temperature was 35 °C. The liquid chromatography (LC) effluent was introduced into the electrospray ionization (ESI) source by a split-flow valve with a ratio of 3:1. All mass spectra were acquired in negative ion mode, and the parameters were as follows: drying gas 4 L/min; drying gas temperature 300 °C; nebulizer (high-purity nitrogen) pressure 15 psi; capillary voltage 4.0 kV; fragmentor voltage 135 V; and cell accelerator voltage 7 V. For full-scan mass spectrometry (MS) analysis, the spectra were recorded in the m/z range of 100–750.

Levels of 29-hydroxy-friedelan-3-one and polpunonic acid were analyzed using Thermo Q Exactive Plus. Chromatographic separation was carried out on a Thermo Syncronis C18 column (2.1 × 100 mm, 1.7 μm). The flow rate of the mobile phase, consisting of 0.1% (v/v) formic acid in water (A) and acetonitrile (B), was set to 0.4 mL/min. The gradient program was as follows: 0–12 min, 10–60% B; 12–17 min, 70% B; 17–25 min, 95% B; 25–28 min, 95% B; 28–29 min, 5% B. Mass spectra were acquired in both positive and negative ion modes with a heated ESI source, and the parameters were as follows: aus. gas flow 10 L/min; aus. gas heater 350 °C; sheath gas flow 40 L/min; spray voltage 3.5 kV; capillary temperature 320 °C. For full-scan MS/data-dependent (ddMS2) analysis, spectra were recorded in the m/z range of 50–750 at a resolution of 17,500 with automatic gain control (AGC) targets of 1 × 106 and 2 × 105, respectively. Levels of metabolites in different tissues were measured by comparing the area of the individual peaks with standard curves obtained from standard compounds.

Yeast in vivo assays were performed following a previously described protocol with some modifications [78]. The yeast expression vector constructs or empty vector were transformed into the yeast Saccharomyces cerevisiae WAT11 [79, 80] using the Yeast Transformation II Kit (ZYMO, CA, USA), and screened on synthetic-dropout (SD) medium lacking uracil (SD-Ura) with 20 g/L glucose. After growing at 28 °C for 48–72 h, transformant colonies were initially grown in 20 ml of SD-Ura liquid medium with 20 g/L glucose at 28 °C for approximately 24 h until the OD600 reached 2–3. Yeast cells were harvested by centrifugation at 4000g and resuspended in 20 mL of SD-Ura liquid medium supplemented with 20 g/L galactose to induce target proteins, while friedelin or 29-hydroxy-friedelane-3-one was applied to the cultures at a final concentration of 25 mM. After 48 h of fermentation (supplemented with 2 mL galactose after 24 h), yeast cells were harvested by centrifugation and extracted with 2 mL of 70% methanol in an ultrasonic water bath for 2 h. The supernatants were filtered with a 0.2-μm Millipore filter and analyzed by LC–MS.

The protocol for enzyme assays in vitro was performed as described previously with some modifications [81]. Yeast transformation and target protein induction were performed as described above, except for 24 h of fermentation after galactose supplementation. Yeast cells were harvested by centrifugation and suspended in a 10-mL mixture of 50 mM Tris-HCl (pH 7.5), 1 mM EDTA, 0.5 mM phenylmethylsulfonyl fluoride, 1 mM dithiothreitol, 0.6 M sorbitol and ddH2O. High pressure cell disruption equipment (Constant Systems, Northants, UK) was used to crush the yeast cells. After centrifugation, approximately 10 mL of supernatant was collected, and CaCl2 was applied at a final concentration of 18 mM. Microsomal proteins were then collected by centrifugation and suspended in storage buffer containing 50 mM Tris-HCl (pH 7.5), 1 mM EDTA and 20% (v/v) glycerol with a final concentration of 10–15 mg/mL determined by the Bradford method [82].

The catalytic activity of putative CYP was assayed in a 100-μl reaction volume, which contained 100 mM sodium phosphate buffer (pH 7.9), 0.5 mM reduced glutathione, 2.5 μg of extracted protein and 100 μM substrate (friedelin or 29-hydroxy-friedelan-3-one). The reaction was initiated by adding NADPH at 1 mM and incubating for 12 h at 28 °C. Methanol was then added at a final concentration of 70% to quench the reaction. The reaction mixture was filtered with a 0.2-μm Millipore filter and analyzed by LC–MS. Microsomal proteins extracted from yeast harboring the empty vector were used as a negative control.

The genomes of T. wilfordii, O. sativa japonica and V. vinifera were compared using MCScan Toolkit v1.1 [83] implemented in Python. The genomes of O. sativa v32 and V. vinifera v32 were downloaded from Ensembl Plants [84]. Syntenic gene pairs were identified using an all-vs-all BLAST search using LAST [85], filtered to remove pairs with scores below 0.7, and clustered into syntenic blocks in MCScan. Microsynteny plots were constructed using MCScan.

We obtained 77.86 Gb of Nanopore reads, amounting to 207.16× coverage of the 375.84-Mb genome, a size estimated by k-mer distribution analysis (Figure 3 and Table 5).

The draft genome was assembled to obtain primary contigs, with a total size of 340.12 Mb and contig N50 of 3.09 Mb (Table 6).

The GC content of the genome was 37.19%, with 0.00% N (Table 7).

Variant calling showed a heterozygosity rate of 0.25% (Table 8).

BUSCO analysis showed 95.2% complete single-copy genes (Table 9).

Short reads obtained from Illumina sequencing in the genome survey were aligned to the genome (Table 5), which exhibited a high consistency with a 95.31% mapping rate and 93.99% coverage (Table 10).

In addition, 87.85% of expressed sequence tags (ESTs) could be identified in the assembly, indicating high coverage of the genome (Table 11).

For assembly improvement, 60.80 Gb (161.77× of estimated genome size) of reads from Bionano sequencing were obtained and integrated with the draft genome to construct scaffolds, which updated the genome of T. wilfordii from a total contig length of 340.12 Mb and a contig N50 of 3.09 Mb to a total scaffold length of 342.59 Mb and a scaffold N50 of 5.43 Mb (Tables 12 and 5).

Furthermore, we anchored 91.02% of the original 342.61-Mb assembly into 23 groups using Hi–C technology (Tables 13 and 14).

All the super-scaffold was able to be placed in one of 23 groups (Figure 4). The super-scaffold N50 reached 13.03 Mb, with the longest super-scaffold being 17.75 Mb in size (Tables 13 and 14). The number of groups, hereafter referred to as pseudochromosomes, corresponded well to the number of chromosomes reported previously [86].

For genome annotation and gene expression profile analyses, roots, stems, young leaves, mature leaves, flower buds and flowers of T. wilfordii plants were collected prior to RNA-seq using the Illumina platform. Furthermore, RNA samples from different tissues were mixed, then sequenced using the PacBio platform to obtain full-length transcriptome sequences (Table 5). A combined strategy involving de novo prediction, homology prediction, RNA-seq and PacBio read alignment was used to construct the gene structure for the T. wilfordii genome. The final set of annotated genes amounted to 31,593 genes, with an average length of 3180 bp and an average coding sequence length of 1182 bp (Table 15).

A total of 27,301 genes (86.41%) were supported by RNA-seq data and 23,229 genes (73.53%) were supported by all the methods used; these genes were annotated with high confidence. Gene function annotation was performed by BLAST analysis of the protein sequences of predicted genes against public databases, including NR, Swiss-Prot, KEGG, GO, Pfam and InterPro. A total of 30,535 (96.70%) gene products could be functionally predicted, and 22,491 sequences could be annotated by at least one of the databases (NR, SwissProt, InterPro and KEGG) (Table 16).

Repeat sequence annotation showed that the T. wilfordii genome contained 44.31% repetitive sequences. Among these sequences, tandem repeats (small satellites and microsatellites) and interspersed repeats accounted for 0.95% and 43.36%, respectively. Long terminal repeats (LTRs) of retroelements were the most abundant interspersed repeats, occupying 36.74% of the genome, including 13.70% Gypsy LTRs and 9.84% Copia LTRs, followed by DNA transposable elements at 1.68% (Table 17).

Noncoding RNA annotation revealed that the T. wilfordii genome possessed 355 microRNAs (miRNAs), 797 transfer RNAs (tRNAs), 827 ribosomal RNAs (rRNAs), and 982 small nuclear RNAs (snRNAs) (Table 18).

Integrated distributions of the genes, repeats, noncoding RNA densities, and all detected segmental duplications are shown in Figure 5.

To identify evolutionary characteristics and gene families, the T. wilfordii genome was compared with 11 published genomes of nine eudicot species (A. thaliana, C. sinensis, V. vinifera, G. max, M. truncatula, G. uralensis, C. sativus, P. trichocarpa and R. communis) and a monocot species (O. sativa). In addition, Amborella trichopoda, one of the basal groups of angiosperms, was selected as an outgroup. Based on gene family clustering analysis, 29,189 gene families were identified, of which 7296 were shared by all 12 species, and 485 of these shared families were single-copy gene families (Figure 6). Gene family numbers were compared between T. wilfordii and four fabid species. As shown in Figure 7A, 10,722 gene families were shared by G. max, C. sativus, G. uralensis, and M. truncatula, and 1086 gene families were specific to T. wilfordii. Compared with the most recent common ancestor (MRCA) of the 12 plant species, 15 gene families with 152 genes were expanded, including CYPs (see GigaDB for table [87]), and 42 gene families with 54 genes were contracted in T. wilfordii (Figure 7B). KEGG analysis showed that the expanded genes were enriched in pathways related to ‘ubiquinone and other terpenoid-quinone biosynthesis’ and ‘steroid hormone biosynthesis’, suggesting that gene family expansion contributed to specialized metabolite biosynthesis in T. wilfordii. A phylogenetic tree was constructed based on the super-alignment matrix of 485 single-copy orthologous genes from the 12 species. The branching order showed that A. thaliana (Brassicales) and C. sinensis (Sapindales) were sister to P. trichocarpa, R. communis (Malpighiales) and T. wilfordii (Celastrales), which diverged approximately 109.1 Mya, followed by divergence of T. wilfordii and species from Malpighiales approximately 102.4 Mya (Figures 5B and 8). These results were consistent with a previously proposed phylogenetic order, in which Celastrales and Malpighiales were found to be sister to each other [88].

The candidate gene identification procedure is illustrated in Figure 9. Based on HMMER analysis, 213 full-length ORFs of CYP genes were extracted from the T. wilfordii genome; these were annotated phylogenetic analysis with CYPs from A. thaliana. Thirty-five CYPs related to triterpenoid oxidases were identified, belonging to different subfamilies, including CYP716, CYP72, CYP71, CYP93, CYP705A and CYP81, which were previously reported to be functionally associated with diverse triterpenoid structural modifications (Figure 10) [89].

On the other hand, the expression patterns of the 213 identified CYPs were identified with TwOSC1 and TwOSC3, which are the two committed enzymes involved in the biosynthesis of the precursor of celastrol in T. wilfordii [21]. Based on RNA-seq data for various tissues, 20 profiles of gene coexpression were obtained, of which only profiles #3 and #13 showed significance (P-value < 0.05) (Figure 11). Profile #3 contained TwOSC3, and 45 CYPs showed similar expression patterns, while profile #13 included TwOSC1, and 51 CYPs had coexpression trends (Figure 12). This suggests that these CYPs are potentially involved in the biosynthesis of celastrol.

To narrow down the candidate genes, the 35 CYPs identified by phylogenetic analysis were compared, and the genes showed patterns of coexpression with TwOSC1 and TwOSC3. As shown in Figure 13A, nine and seven triterpenoid biosynthesis-related CYPs showed patterns of coexpression with TwOSC1 and TwOSC3, respectively. However, no CYPs were common between the TwOSC1 group and TwOSC3 group. Based on tissue expression profiles, the 16 CYPs were clustered separately into two clades with TwOSC1 and TwOSC3, in which TwOSC3 exhibited root-specific expression, while TwOSC1 was highly expressed in leaves and other aerial parts (Figure 13B).

Gene-to-gene and gene-to-metabolite Pearson’s correlation coefficients (r) were calculated using the tissue expression profiles of the 16 outstanding CYPs mentioned above, as well as three other known genes related to celastrol biosynthesis (TwHMGR1, TwFPS1 and TwDXR) [90–92], and the known intermediate product and celastrol concentrations [21]. As shown in Figure 13C, seven CYPs positively correlated with celastrol biosynthesis-related genes, with high Pearson’s r and significant P-values. In addition, these CYPs highly correlated with the levels of 29-hydroxy-friedelan-3-one, polpunonic acid, wilforic acid A and celastrol, which all specifically accumulated in the roots of T. wilfordii (Figure 14).

Phylogenetic analysis placed these seven CYPs and CYPs from other species into three clades representing different functions in the structural modifications of triterpenoids (Figure 13D). Two CYPs (tw18g03520.1 and tw01g08960.1) were clustered with CYP81Q58 from C. sativus, which catalyzes hydroxylation of the C-25 position in cucurbitacins [93]; four CYPs (tw12g11350.1, tw18g06010.1, tw18g06210.1 and tw23g06640.1) were clustered with CYP712K4 from Monteverdia ilicifolia, which catalyzes the oxidation of the C-29 position using friedelin as a substrate [94]; and tw03g08450.1 was clustered with CYP716C11 from Centella asiatica, which hydroxylates the C-2 position of oleanolic acid and ursolic acid [95].

Since we were interested in identifying a C-29 position oxidase that could catalyze the conversion of friedelin to polpunonic acid, 3 CYPs were finally chosen as candidates for functional validation: tw18g06010.1, tw18g06210.1 and tw23g06640.1 (hereafter referred to as TwCYP712K1, TwCYP712K2 and TwCYP712K3) according to the closer relationship with CYP712K4 from M. ilicifolia, which is related to C-29 hydroxylation [94].

The full-length ORFs of TwCYP712K1, TwCYP712K2 and TwCYP712K3 were successfully clones and separately expressed in yeast fed with friedelin or 29-hydroxy-friedelan-3-one. However, no new peaks could be detected from the yeast strains expressing the enzymes and supplemented with friedelin compared with the empty vector (EV) control (see GigaDB [87]). This probably because the hydrophobic substrate could not be transported into the yeast cells. When fed with 29-hydroxy-friedelan-3-one, both TwCYP712K1 and TwCYP712K2 converted the substrate to a new compound possessing the same mass charge ratio (m/z) as the polpunonic acid standard, while TwCYP712K3 showed no such activity in this assay (Figure 15A and C). To further explore the enzyme activities, microsomes were extracted from yeast cells and the proteins incubated with friedelin or 29-hydroxy-friedelan-3-one for 12 h. As shown in Figure 15B, TwCYP712K1 and TwCYP712K2 converted 29-hydroxy-friedelan-3-one to a new compound in vitro, consistent with previous yeast in vivo assays. In addition, no new peak was detected from the enzyme reactions supplemented with friedelin compared with the EV control (see GigaDB [87]).

Genome analysis showed that TwCYP712K1 (tw18g06010.1) and TwCYP712K2 (tw18g06210.1) were located in pseudochromosome 18 within an approximately 200-kb region. This led us to examine the evolutionary relationship between specific CYPs. We hypothesized that a gene duplication event must have occurred during evolution. However, amino acid alignment indicated only 70.57% identity between TwCYP712K1 and TwCYP712K2 (Figure 16), suggesting that these two genes became specialized a long time ago. Syntenic analysis showed that TwCYP712K1 and TwCYP712K2 had corresponding collinear genes in V. vinifera but not in O. sativa, indicating that these two genes appeared after the species differentiation of Poaceae and Vitales (<178.6 Mya), and they came from the common ancestor but divided after the species differentiation of Vitales (<130.6 Mya) (Figures 7 and 17).