E. pubescens plants were grown in the germplasm nursery at the Shawan District (29°N, 103°E), Leshan city, Sichuan province of China. In the germplasm nurseries, Epimedium plants were grown and covered with a black shade net to avoid direct sun exposure. Wild populations of E. pubescens were thoroughly investigated during field trips from 2018 to 2019 by Dr. Chaoqun Xu ( Table S2 ) and totally 39 samples of E. pubescens were identified by Professor Baolin Guo. Voucher specimens from these wild populations were deposited at the herbarium of the Institute of Medicinal Plant Development (IMPLAD), Beijing, China, under the voucher numbers from B. L. Guo00841 to B. L. Guo00879. Fully developed leaves were detached in the spring of 2018 and washed with double-distilled sterile water, flash frozen in liquid nitrogen and stored at -80°C until DNA extraction and sequencing. Fresh roots, shoots, leaves, flowers, and fruits from the same individual plant were harvested and immediately frozen in liquid nitrogen, and stored at -80°C for RNA extraction and RNA sequencing.

Genomic DNA was extracted from the leaves using the CTAB DNA extraction protocol (Varma et al., 2007). DNA concentration and purity were examined using NanoDrop and Qubit (Thermo Fisher Scientific, MA, USA), and DNA integrity was assessed using the pulsed-field electrophoresis. The genomic DNA was then used to construct libraries with an average insert size of 300 bp, and sequenced on the Illumina Novaseq 6000 platform (Illumina, San Diego, CA, USA) with average sequencing depth at about 230.94 × (771 Gb). For ONT sequencing, the high-quality genomic DNA was separated using the BluePippin™ System (Sage Science, USA), and purified to construct a library with size ranging from 15 to 50  kb using ONT template prep kit (SQK-LSK109, Oxford headquarters, USA) and NEB Next FFPE DNA Repair Mix kit (New England BioLabs, MA, USA). The high-quality library was sequenced on the ONT PromethION platform with MinION flow cell (R9.4.1) and ONT sequencing kit (Oxford headquarters, USA. Finally, a total of 906 Gb raw data (271.32 ×) was generated with an average ONT length of 19.79 kb and an N50 of 31.17 kb ( Table S3 , S4 ). Raw data was processed for base calling by Oxford Nanopore base caller using default parameters (Wick et al., 2019), and ONT reads were trimmed with Porechop (https://github.com/xxz19900/Porechop). Finally, these reads were assembled into contigs using SMARTdenovo (Liu et al., 2021a). The assembled contigs were then polished with Illumina short reads three times by Racon (Vaser et al., 2017) and Pilon v1.20 (Walker et al., 2014). Finally, BWA-MEM2 (https://github.com/bwa-mem2/bwa-mem2) for alignment of short reads and BUSCO 4.0.6 with viridiplantae_odb10 (Simao et al., 2015) were used to assess the quality and completeness of the assembly.

For the prediction of gene models, ab initio prediction program Augustus v2.4 and evidence/homology-based strategies were applied to annotate genomic contigs and transcriptomic data and the results were integrated into final gene models using EVM v1.1.1 (Haas et al., 2008). For repeat annotation, the sequences of genome assembly were subjected to structural and ab initio prediction of repeats using LTR_FINDER v4.0.6, RepeatScout v1.0.5 and RepeatMasker v4.0.6. For non-coding RNA annotation, microRNA, rRNA, tRNA, and other functional RNA were predicted by combining several strategies. For functional annotation, the predicted gene models were subjected to homology searches against the following databases: NCBI non-redundant protein sequences (NR), Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups of proteins (KOG/COG/eggNOG), Gene Ontology (GO), Pfam and Swiss-Prot/TrEMBL ( Table S5 ). In the RNA-seq analysis, raw data was cleaned by Trimmomatic software v0.39 and the gene expression levels were determined with Hisat2 and Stringtie. The differentially expressed genes were detected by R package DESeq v1.10.1. Co-expressed gene network was inferred from the FPKM values in different tissues (FPKM > 10) using weighted gene co-expression network analysis (WGCNA) (Langfelder and Horvath, 2008).

Protein sequences from E. pubescens and other 11 angiosperms ( Table S6 ) were clustered into orthologous groups using Orthofinder software v2.3.3 and aligned by MUSCLE v3.8.155147. The phylogenetic tree was built using the maximum-likelihood method with 1,000 bootstrap replicates in RAxML v8.2.1248. Synonymous substitution rate per synonymous site (Ks) was calculated by MCScanX software. Ks peaks were determined by using the Genome_tools (Ks_Density_plot.r, https://github.com/ZhangXu-CAS/Genome_tools/), and the Ks distribution plots were made using R package ggplot2. The times of Ks peaks were calculated by the formula, Ks/13×1000 (MYA) (Gaut et al., 1996). The divergence times among different taxonomic groups were estimated by the TimeTree online service (Kumar et al., 2017).

To identify the candidate prenyltransferase genes, protein sequences from the known plant flavonoid prenyltransferases were used to search the E. pubescens genome using BLASTP with E-value of 10⁻⁵. The candidate sequences were further submitted to NCBI CDD database to confirm the presence of conservative domains. The protein properties of putative EpPTs, including physical and chemical properties, subcellular localization, transmembrane (TM) α-helices and the presence of chloroplast transit peptides (cTP) were predicted using the online ExPASY ProtParam tool (http://web.expasy.org/protparam/), TargetP online server (http://www.cbs.dtu.dk/services/TargetP/), TMHMM 2.0 (http://www.cbs.dtu.dk/services/TMHMM/), ChloroP 1.1 (http://www.cbs.dtu.dk/services/ChloroP/) and PSORT (http://psort1.hgc.jp/form.html), respectively. The sequence alignments were generated using ClustalW (Li, 2003), and the gene tree was constructed using RAxML package v 8.13 with 1,000 bootstrap replicates (Stamatakis, 2014). The best substitution model was determined using BestModel.

Total RNA from fresh leaves of E. pubescens was extracted using an Eastep^® Super total RNA Extraction Kit (Promega, Shanghai, China). First-strand cDNAs were synthesized using FastKing One-Step RT-PCR Kit (TIANGEN Biotech, Beijing, China) for amplification of putative flavonoid E. pubescens prenyltransferase genes (EpPTs). Specific primers ( Table S7 ) were designed according to the candidate gene sequences and nested PCR amplification was performed using Q5^® High-Fidelity DNA Polymerases (New England BioLabs, MA, USA). The PCR program was set up as follows: denaturation at 98°C for 30 s; 35 cycles of 98°C for 10 s, 53-56°C for 30 s, 72°C for 1 min and 20 s; and a final extension at 72°C for 5 min. After the first-round PCR amplification, PCR products were used as the template for second-round PCR amplification under the same PCR conditions. PCR products were purified using AxyPrep DNA Gel Extraction Kit (Corning, NY, USA) and further cloned into pTOPO-Blunt simple vectors (LANY, Beijing, China), which were transformed into Escherichia coli competent cell Trans1-T1 (TransGen Biotech, Beijing, China) and confirmed by Sanger sequencing. Candidate EpPT genes and two truncated EpPTΔTP constructs were cloned into the entry vector pENTR/D-TOPO (Invitrogen, CA, USA) and confirmed by Sanger sequencing, and then they were inserted into yeast expression vectors (pDR196GW) using LR Clonase™ II Enzyme (Invitrogen, CA, USA). The resulting vectors, pDR196GW-EpPT and pDR196GW-EpPTΔTP were separately transformed into yeast strain DD104 using the modified LiAc method (Pompon et al., 1996; Liu et al., 2003). The transformants were screened on SD/-Ura plates and confirmed by PCR. Yeast expression vector pDR196GW and strain DD104 were kindly provided by Professor Guodong Wang (Institute of Genetics and Developmental Biology, the Chinese Academy of Sciences, Beijing, China). Yeast-Extract Peptone Adenine Dextrose Medium (YEPAD) and Ura Minus Medium were purchased from FunGenome (Beijing FunGenome Co. Ltd, Beijing, China).

Monoclonal yeast clones, which contained empty pDR196GW vector (blank control) or pDR196GW-EpPTs/EpPTΔTPs vector, were cultured in 5 mL of SD (-Ura) overnight at 28°C, separately. Cultured yeast of 200 μL was inoculated into 780 μL of SD (-Ura) broth supplemented with 17 flavonoid substrates ( Figure S3 ) with a final concentration of 200 μM and grown at 28°C for 72 h. These flavonoid substrates were purchased from Shanghai Yuanye Bio-Technology Co. Ltd (Shanghai, China). After the incubation period, the enzymatic reaction mix was ultrasonically extracted three times with an equal volume of ethyl acetate for 20 min. The ethyl acetate solvent in samples was evaporated and the dried powder was then dissolved in 200 μL absolute methanol for HPLC and UHPLC-PDA-Q-TOF/MS analyses.

The yeast cultures fed with kaempferol were scaled up to 1.5 L, and extracted with ethyl acetate. For the isolation of enzymatic products, semipreparative RP-HPLC was conducted on a Lumtech K-501 equipped with a YMCPack ODS-A column (250 mm × 10 mm i.d., 5 mm, YMC Co., Ltd., Kyoto, Japan) at a flow rate of 3 mL·min^-1. The solvent system consisted of a linear gradient (70%–100%, v/v) of methanol in water over 0-20 min. UV detection was set at 254 nm and 280 nm.

LC-MS/MS analysis was performed on Waters Xevo G2-XS Tof (Waters, Milford, MA, USA). The separation was carried out with a Waters ACQUITYTM HSS T3 C18 column (2.1 mm×100 mm, 1.8 μm) at 40°C. The gradient is consisted of 0.1% formic acid (A) and acetonitrile (B) as the mobile phase, 0-1.5 min (21%-24% B), 1.5-3 min (24%-25% B), 3-4 min (25%-29% B); 4-5 min (29% B); 5-6.5 min (29%-32% B); 6.5-7 min (32%-44% B); 7-8 min (44%-45% B); 8-9 min (45%-46% B); 9-11min (46%-95% B). The operating conditions were as follows: flow rate of 0.6 mL·min^-1 with positive ion ESI mode, a capillary voltage at 3 kV, a cone voltage at 50 V, desolvation gas with a flow rate of 850 L·h^-1. The mass-to-charge ratio was scanned from 100 to 1,600 m/z.

Approximately 10 mg of each compound was evaporated to dryness under N₂ gas, resuspended in dimethyl sulfoxide-d ₆ (DMSO-d ₆), and analyzed through ¹H NMR spectra acquired on a Bruker 600 spectrometer (Bruker, Rheinstetten, Germany).

¹H NMR data of substrate and prepared isoprene products were shown as follows:

Both previous studies and our survey showed that the genome size of Epimedium species varied from 3.14 Gb/1C to 4.49 Gb/1C, and the size of E. pubescens genome was estimated to be 3.23 Gb based on k-mer (k=23) distribution analysis and flow cytometry analysis (Chen et al., 2012; Liu et al., 2013; Zhang et al., 2018) ( Table S1 ). In addition, the k-mer analysis of E. pubescens further revealed a relatively low hybridization rate at 1.2% and an estimated percentage of repetitive elements at 61.57%, suggesting that E. pubescens possessed a relatively noncomplicated genome (the E. pubescens genome is complicated compared to model diploid plants such as Arabidopsis thaliana and Oryza sativa, but uncomplicated compared to allopolyploid plants Triticum aestivum and Medicago sativa) (Michael, 2014). Therefore, E. pubescens was eventually selected as a representative species for the construction of reference genome for Epimedium genus.

Approximately 48.01 million Oxford Nanopore Technology (ONT) long reads were acquired and accounted for 906 Gb (271.32×) with an average read size of 20,267 bp ( Table S3 , S4 ). The preliminary genome assembly of E. pubescens was created with 39,251,621 ONT reads of >2 kb in length, yielding 6,229 contigs with a total size of 3.34 Gb ( Table S8 ). Illumina reads of 771 Gb (230.94×) was used for further polishing the initial genome assembly to achieve contig N50 of 871.79 kb ( Table 1 and Table S8 ). Subsequently, 1,965,567,758 paired-end reads from the Hi-C (High-throughput chromosome conformation capture) sequencing were used to successfully anchor 6,176 contigs (99.15%) onto six chromosome-level pseudo-molecules, with an average size of 538.98 Mb ( Figure 2A and Table S9 ). The resultant chromosome-level genome assembly had a size of 3.34 Gb with evenly distributed genes, slightly unevenly distributed Transposable Elements (TEs) and 37.48% of average GC content ( Table 1 and Figure 2B ). Eventually, 5,162 contigs of 3.16 Gb were fully anchored and oriented on six complete chromosome-level pseudo-molecules ( Table 1 and Figure 2B ).

Homology and ab initio based gene prediction strategies were combined to predict 44,722 protein-coding gene models from E. pubescens genome assembly ( Table 1 ). Among these gene models, protein coding genes accounted for 487.81 Mb (14.61%) with an average gene length of 13.36 kb ( Table 1 ), and have 33,355 with 1~5 exons (74.58%) and 11,367 with more than 5 exons (25.40%) ( Table S10 ). A total of 94.5% gene models were successfully annotated by eight protein databases ( Table S5 ), and a total of 8,415 noncoding RNA genes were identified, including 1,713 rRNA (ribosomal RNA), 172 miRNA (microRNAs), 5,794 tRNA (transfer RNA), 278 snoRNA (small nucleolar RNA), and 458 snRNA (small nuclear RNA) ( Table S11 ). In addition, large amounts of repetitive elements were identified in the E. pubescens genome (66.93%), which was relatively high in comparison with other species of Ranunculales ( Table S12 ). Among these repetitive elements, there were 62.61% Class I retrotransposons, which were predominantly the long terminal repeats (LTRs) with mainly Gypsy and Copia at 26.62% and 5.03% respectively, and 7.92% of DNA transposons (Class II) ( Table S12 , S13 ). Notably, 1,236 out of 1,375 (89.9%) BUSCO (Benchmarking universal single-copy orthologues) core genes were confirmed in the final chromosome-level E. pubescens genome ( Table S14 ). These above evidences indicated that the current chromosome-level genome assembly of E. pubescens was relatively complete and accurate.

The genomic syntenic blocks and orthologous gene ratio between V. vinifera and E. pubescens was analyzed, demonstrating that there was a significant percentage of genes with 3:1 ratio of V. vinifera to E. pubescens, but not be observed between P. somniferum/C. chinensis and E. pubescens ( Figures 3B, C and Figure S1A–C ). It was known that only one Whole Genome Triplication (WGT-γ) event occurred in Vitis vinifera (Jaillon et al., 2007), and moreover, P. somniferum and C. chinensis were also proven to have not experience WGT-γ event (Guo et al., 2018; Liu et al., 2021c). The combination of previous research results with current evidences suggested that E. pubescens had escaped the ancient WGT-γ event. To further dissect genome evolutionary history, the analysis of synonymous substitution rate per synonymous site (Ks) was performed among the orthologous genes of E. pubescens, P. somniferum, Coptis chinensis and Aquilegia coerulea from Ranunculales, revealing an significant peak (around Ks=0.997) in E. pubescens compared to C. chinensis (around Ks=0.85) ( Figure 3A ), indicating that only one Whole Genome Duplication (WGD) event occurred during E. pubescens genome evolution, which was consistent with the widespread occurrence of WGD events in flowering plants (Jaillon et al., 2007).

To investigate the evolving history of E. pubescens, a phylogenetic tree was constructed using shared single-copy orthologues identified by Orthofinder from E. pubescens and 11 key angiosperm species, including A. trichopoda, P. somniferum, C. chinensis and V. vinifera, etc. ( Figure 3D and Table S6 ). In the resultant phylogenetic tree, E. pubescens and other species from Ranunculales formed an early diverging taxonomic clade of basal eudicots, and these species were further divided into 3 distinct taxonomic groups, including Ranunculaceae (A. coerulea and C. chinensis), Papaveraceae (P. somniferum and Macleaya cordata) and Berberidaceae (E. pubescens). Using Treetime online services, the divergence time between Ranunculaceae and Berberidaceae was determined at ~81 MYA (million years ago) and the formation of Ranunculaceae was also estimated at ~66 MYA. In addition, there exhibited some strong patterns of inter-chromosomal synteny between chromosomes of E. pubescens, as shown in Figure S1C , inferring that multiple ancient chromosomal breaks and fusion events could occur after the only WGD event in E. pubescens genome.

It has been known that, for Epimedium flavonoids, after the formation of basic flavonoid skeleton, several enzymatic steps still are needed for further modification, including prenylation, methylation and glycosylation (Pandey et al., 2016; Nabavi et al., 2020). The addition of prenyl group to flavonoids is one of the most important enzymatic steps, which is catalyzed by prenyltransferases (PTs) ( Figure 4A ), and provides the medicinal efficacy for Epimedium flavonoids (Ming et al., 2013; Mbachu et al., 2020). Genome-wide homology-based search discovered 19 potential PTs from UbiA (ubiquinone biosynthesis gene A) superfamily in E. pubescens ( Figure 4B and Table 2 ). Subsequently, a phylogenetic tree was constructed with these 19 EpPTs and additional 108 UbiA PTs from other plant species ( Table 2 and Table S15 ), revealing that all PTs of UbiA superfamily could be clustered into distinctive groups based on the preference of their substrates, such as flavonoids, polyphenol, chlorophyllide a/b and homogentisate acid etc. ( Figure 4B ). It was found that, 11 EpPTs formed a distinctive cluster that was nested in the clade of PTs using polyphenols or flavonoids as preferred substrates, and six out of these 11 EpPTs (EpPT3, EpPT4, EpPT5, EpPT6, EpPT7, EpPT8 and EpPT9) were found to form a cluster at the end of Chr02 (data not shown).

To further explore these 11 E. pubescens PTs, their tissue expression profiles were analyzed, showing that EpPT8 was most highly expressed in leaves, where the bioactive compounds accumulated ( Figure 4C ) (Guo and Xiao, 1996; Zhou et al., 2012). Further analysis revealed that EpPT8 shared relatively high sequence similarity with known plant flavonoid PT genes ( Figure S2 ), suggesting that EpPT8 was a promising candidate gene for the addition of prenyl group on flavonoids in E. pubescens. In EpPT8, the transit peptides were predicted to be 32 or 82 amino acids in length at the N terminus which were subsequently removed to generate the truncated EpPT8 constructs. Yeast strain DD104 cells were transformed with EpPT8 and two truncated EpPT8 constructs, and subsequently co-cultured with 17 representative flavonoid substrates listed in Figure S3 . After purification, the reaction mixtures were subjected to the examination of PT enzymatic activities. Notably, the recombinant EpPT8 proteins only showed enzymatic activity towards kaempferol (highest enzymatic activity), apigenin and quercetin out of 17 substrates ( Figure S3 ). Meanwhile, the other two truncated EpPT8 proteins have similar levels of prenylation activities ( Table 3 ), suggesting that putative transit peptides in EpPT8 did not affect the enzymatic activity of EpPT8 proteins.

To examine the products in the EpPT8 catalyzed reaction using kaempferol as substrate, their HPLC chromatographs, MS spectra and MS/MS spectra were acquired for the confirmation of chemical property. The above reaction products were shown at about 55% conversion ratio of substrate to product and possessed the same HPLC retention time and reference data as 8-prenylkaempferol standard ( Figure 4D , Figure S3 , S4 ) (Liu et al., 2021b), indicating that EpPT8 was able to prenylate kaempferol at C-8 position. To confirm the chemical structure, the EpPT8 catalyzed reaction products were further separated with preparative liquid chromatograph and determined by NMR experiments. In ¹H NMR spectrum ( Figure S5 ), the characteristic signals (Hillerns and Wink, 2005; Kim et al., 2018) were identified as follows: δ _H=1.63 ppm (3H, s, H-5′′) and δ _H=1.75 ppm (3H, s, H-4′′) for two methyl groups, δ _H=3.43 ppm (2H, d, J=7.0 Hz, H-1′′) for one methylene and δ _H=5.18 (1H, t, J=6.9 Hz, H-2′′) for one methine, suggesting the addition of a dimethylallyl moiety onto kaempferol (De Souza et al., 2017). By comparing the ¹H NMR spectrum of kaempferol (De Souza et al., 2017) with that of 8-prenylkaempferol in the literature (Hillerns and Wink, 2005), the addition of dimethylallyl unit was determined to occur at the C-8 position ( Figure S5 ). Taken together, EpPT8 could catalyze the prenylation of kaempferol, leading to the principal enzymatic product, 8-prenylkaempferol ( Figure 4D ). Moreover, the enzymatic products of EpPT8 with quercetin and apigenin as substrates were also predicted to be prenylated products with expected MS and MS/MS spectrum ( Figures 5A, B ), but the prenylation position could not be resolved by NMR due to an extremely low amount of available enzymatic products. In addition, among 11 PTs, EpPT1, EpPT4 and EpPT9 were close to full length with minor change (similarity 88.52%-95.4%) compared to EpPT8, but they are either failed to be cloned or proved no PT activity. The rest of the PTs are with low similarity to EpPT8 and either N-terminal or C-terminal truncated.