A chromosome-level Pinellia ternata genome assembly provides insight into the evolutionary origin of ephedrine and acrid raphide formation

Pinelliae Rhizoma, the dried tuber of Pinellia ternata (Thunb.) Breit. (Araceae), is a famous bulk materia medica broadly used in traditional herbal medicines in Asian countries^[1−3]. Particularly, Xiaoqinglong Tang (Shoseiryuto), Banxia Houpu Tang (Hangekobokuto), and Wenjing Tang (Unkeito) are the classical formulas widely used in Traditional Chinese Medicine (TCM) and Japanese Kampo medicine; the dried tuber of P. ternata is the key component of these medicines^[4]. Moreover, other formulas, such as Banxia Shumi Tang, are also broadly applied in TCM. Medicines made using P. ternata primarily have antiemetic, expectorant, sedative, antipyretic, and styptic effects^[5−8]. These medicinal functions are attributed to various bioactive compounds, including alkaloids, iridoid glycosides, and organic acids in the tuber, as reported in several pharmacological studies^[9−12]. Of the compositions, alkaloids are deemed as the core component of the medicinal efficacy of P. ternata tubers^[12−14].

The alkaloids present in P. ternata are guanosine, inosine, trigonelline, and ephedrine^[15]. Among them, ephedrine is of particular interest due to its function as an antitussive^[16,17]. Ephedrine, the key bioactive component of P. ternata, accumulates predominantly in tubers^[18], and the ephedrine content in tubers typically ranges from 6 to 14 μg·g⁻¹ dry weight^[19]. In addition to expanding the cultivation acreage, the breeding or selection of P. ternata with high medicinal components is an important approach to meet the market demand. Understanding the biosynthesis pathway of ephedrine may allow for an increase in its accumulation and a decrease in the dose of P. ternata tuber. Previous studies have investigated the biosynthesis pathway of ephedrine in P. ternata. Specifically, Zhang et al.^[20] and Duan et al.^[14] conducted a transcriptome analysis and external elicitation to construct the biosynthetic pathway of ephedrine in P. ternata. However, our understanding of the pathway was based on transcriptome data, which is insufficient to provide complete candidate genes for each biosynthatase. In addition, the final steps from 1-phenylpropane-1,2-dione to ephedrine remain unclear.

In addition to its medically relevant alkaloids, P. ternata has an acrid taste and is defined as 'toxic' in The Chinese Pharmacopoeia^[21]. The acridity of P. ternata tubers, as well as the acridity of several other Araceae species is mainly attributed to insoluble needle-like crystals called acrid raphides^[22,23]. Recent studies showed that the raphides are composed of calcium oxalate, proteins, and polysaccharides^[24]. Acrid raphides are calcium oxalate by a conjugated reaction between calcium and oxalic acid, and exogenous fertigation of calcium could enhance calcium oxalate^[25]. The acridity of the P. ternata tuber can be explained by the calcium oxalate crystals that mediate the introduction of inflammatory proteins (lectins) by penetrating the mucous membrane^[26]. Because of the presence of acrid raphides, proper processing via decocting with licorice, alum, or ginger is required while making traditional medicines using P. ternata tubers^[4,27]. Processing, however, would lower alkaloid concentrations significantly in P. ternata tubers^[28]. Though the composition of raphides is understood, it is currently unclear how acrid raphides form in Araceae species.

Though there are approximately 110 genera and 3,500 species in Araceae, only Spirodela polyrhiza^[29], Colocasia esculenta^[30], Pinellia pedatisecta^[31], Pistia stratiotes^[32] and Amorphophallus konjac K. Koch^[33] have fully sequenced genomes. Moreover, S. polyrhiza and P. stratiotes are an aquatic plant, and C. esculenta and A. konjac are vegetable crops; therefore, their genome sequences may not be sufficiently compatible to provide a strong reference for P. ternata. P. pedatisecta is a widely used as traditional medicinal herbs^[9]. However, P. ternata and P. pedatisecta belong to different species, which will be different in many aspects. The lack of genome information largely restricts molecular studies on P. ternata. Typically, P. ternata is a hexaploid or octoploid in production^[34], increasing the complexity of genome sequencing. However, using a wild diploid accession, we sequenced the P. ternata genome and provided a chromosome-scale assembly by combining Illumina, Pacbio HiFi, and Hi-C sequencing technology. Further, the biosynthesis pathways for amphetamine-type alkaloids and acrid raphides were constructed. The data presented herein will provide valuable information for illustrating the molecular mechanisms controlling the biosynthesis of alkaloids and acrid raphides.

Individual P. ternata Breit. plants that were collected from their wild habitat in Baojia Town, Yuexi County, Anhui Province, PR China (30°83′ N, 116°08′ E) and identified as diploids via flow cytometry, were planted in the greenhouse at Huaibei Normal University. High-quality genomic DNA was extracted from young leaves of P. ternata for genome sequencing. A paired end library was constructed and sequenced using the BioMarker Technology Company (BioMarker, Qingdao, China) and Illumina NovaSeq X platforms to generate short reads. The PacBio sequel II platform was used for Pacbio HiFi library constructing and sequencing, thus obtaining long reads. The Hi-C library was constructed using the same leaf samples as described previously and sequenced using Illumina NovaSeq X platforms. Hi-C-Pro was applied for quality control and the identification of valid interaction pairs.

Genome size estimation

Two methods were used to estimate the genome size of P. ternata: flow cytometry and K-mer counting. GenomeScope2.0^[35] was used to estimate the characteristics of the P. ternata genome, including genome size, heterozygosity, and repetitiveness based on the K-mer method.

Genome assembly and quality assessment

Pacbio HiFi long reads were assembled using HiFiasm^[36]. Redundant sequences in the initial assembly of P. ternata were removed using Khaper (https://github.com/lardo/khaper). NextPolish^[37] was used to polish the genome assembly, with Illumina reads employed to generate the final contig-level assembly. ALL Hi-C was used to anchor the contigs to chromosomes based on the HiC reads, and JUICEBOX was used to correct assembly errors^[38], allowing for the attainment of a chromosome-level genome. BUSCO (v5.2.2)^[39] with OrthoDB (embryophyta_odb10), Illumina short reads and transcriptome data were mapped to the P. ternata genome with BWA-mem to assess the quality of the assembled genome. The Long-Terminal Repeat (LTR) Assembly Index (LAI)^[40] was also used to evaluate the quality of the genome. The evaluation criteria were as follows: the genome assembly will be considered a draft genome if the LAI value was between 0 and 10, the genome assembly will be considered a reference genome if the LAI value was between 10 and 20, and the genome assembly will be taken as gold if the LAI value was greater than 20.

Gene prediction and functional annotation

GETA, a homology-based tool integrating three methods, was used to predict the protein-coding genes in the P. ternata genome (https://github.com/chenlianfu/geta). Genes were functionally annotated by searching against the eggNOG (Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups) 5.0 database. The protein-coding sequences were submitted to the EggNOGmapper software online with the default parameter settings.

Repetitive sequence annotation

RepeatMasker (v4.0.9) and RepeatModeler (v1.0.8)^[41] were used to identify repetitive sequences in the P. ternata genome. RepeatModeler was used for de novo prediction and a library of consistent repetitive sequences was established. RepeatMasker was then used for comprehensive annotation by loading the repetitive sequence library. Unknown repeat sequences were further classified using TEclass^[42]. Long-terminal repeat retrotransposons (LTR-RTs) were identified using LTR Finder (v1.0.2)^[43] and LTRharvest.

Genome evolution and divergence time estimation

To conduct an evolutionary analysis, we collected the protein sequences of P. ternata and 13 other species: Arabidopsis thaliana, Carica papaya, Oryza sativa, Brachypodium distachyon, Areca catechu, Wolffia australiana, Elaeis guineensis, Cocos nucifera, Colocasia esculenta, Amorphophallus konjac, Pinellia pedatisecta, Pistia stratiotes, and Zostera marina. We also identified the gene clusters of these genomes using Orthofinder^[44] with default settings. CAFE (v4.2.1)^[45] was used to identify the gene families that underwent expansion or contraction in the 11 sequenced species.

Single-copy orthologous genes from the 14 species were extracted and the protein sequences were aligned using MAFFT (v7.307)^[46]. RaxML (v 8.2.12)^[47] was used to construct the maximum-likelihood phylogenetic tree. The MCMCTree program in the PAML package^[48] was used to estimate the species divergence times using the divergence time between A. thaliana and C. papaya (i.e., 68–72 million years ago) and the monocot and eudicot divergence time (i.e., 120–140 million years ago), as shown by A. thaliana and O. sativa as calibrators.

WGD event andsSynteny analysis

P. ternata, A. konjac, P. stratiotes, and P. pedatisecta belong to the subfamily of Aroideae of the Araceae family. To identify potential WGD events of the P. ternata genome, the WGD pipeline^[49] was used to calculate the distribution of synonymous substitutions per synonymous site (Ks).

JCVI^[50] with default parameters was used to identify collinear blocks between P. ternata, A. konjac, and P. pedatisecta, and a dot plot was drawn to confirm synteny between P. ternata and P. pedatisecta genomes.

Integrated genomic and transcriptomic analysis of ephedrine and oxalate biosynthesis-associated genes

Total RNA was extracted from tubers during different growth periods (young (Y), medium (M), and old (O)), as well as roots, stems, leaves, and flowers using the Plant RNA Kit following the manufacturer's instructions. Three biological replicates were conducted for each tuber sample. The isolated RNA was assessed using a NanoDrop Qubit 2.0 Fluorometer. High-quality RNA was sent to Berry Genomics (www.berrygenomics.com) for library construction and sequencing.

The sequencing reads were aligned to the assembled genome using HISAT2^[51]. Stringtie2^[52] was used to calculate the expression levels of all genes through fragments per kilobase of transcript per million mapped fragments (FPKM). DeSeq2^[53] was used to analyze the significantly differentially expressed genes (DEGs) with a false discovery rate (FDR) of ≤ 0.05 and an absolute log2 (fold change) of ≥ 1 as the threshold. DEGs were clustered into different expression files using R package ClusterGVis^[54]. iTAK (v 1.7)^[55] was used to identify the transcription factors (TFs) in the genome and classify them based on the PlnTFDB and PlantTFDB databases.

After removing redundant sequences, we identified genes potentially involved in the biosynthesis pathway of ephedrine and oxalate using BLASTP (E-value =1e⁻⁵)^[56].

According to the flow cytometry results, the estimated genome size of the P. ternata genome (2n = 2x = 26) was ~1.9 Gb (Supplemental Fig. S1, Supplemental Table S1). Based on K-mer counting, a similar genome size was obtained (Supplemental Fig. S2, Supplemental Table S2). The heterozygosity and repetitiveness of the P. ternata genome were 2.51% and 76.2% (Supplemental Table S2) respectively, indicating that the P. ternata genome was a complex diploid genome with high heterozygosity and high repetition. A total of 101.61 Gb Illumina short reads (~48.8×), 65.78 Gb Pacbio HiFi long reads (~31.6×), and 216.46 Gb Hi-C reads (~104.1×) were obtained for the assembly of the P. ternata genome (Supplemental Tables S3−S5). The final assembled sequence was 2.08 Gb with anchoring on 13 pseudo-chromosomes that ranged from 117.77 Mb to 189.55 Mb long (Fig. 1, Supplemental Table S6 & S7). The total length of the 13 pseudo-chromosomes accounted for 95.46% of the total P. ternata genome sequence (Supplemental Fig. S3, Supplemental Table S7). The BUSCO score was 92.60% (1,495), and 84.8% (1,369) of the genes were considered to be complete and single copy (Supplemental Table S8). The LAI value was 21.87 (LAI > 20) (Supplemental Table S9).

Figure 1. 

Genome features of Pinellia ternata. (a) Chromosomes (chr01-Mbchr13); (b) Gene density; (c) DNA transposable elements; (d) LTR/Copia; (e) LTR/Gypsy; (f) P. ternata.

We predicted 34,342 protein-coding genes in the P. ternata genome by integrating homology-based, de novo, and transcriptome-based methods. Among the predicted genes, 32,140 (93.6%) genes were functionally annotated using EggNOGmapper. BUSCO (v5.2.2) analysis found that 1,332 (82.5%) were identified as complete (Supplemental Table S8).

Through a combination of de novo and homology-based approaches, 86.02% of the P. ternata genome assembly was annotated as repeat sequences, including 72.36% of retrotransposons and 10.9% transposons (Supplemental Table S10). The LTR retrotransposons accounted for the largest proportion (62.59%) of retrotransposons and their length was approximately 1,243.8 Mb, including 342.08 Mb of Ty1/Copia, 722.48 Mb of Ty3/Gyspy, and 179.24 Mb of other retrotransposons, which accounted for 17.21%, 36.36%, and 9.02% of genome assembly (Supplemental Table S10), respectively.

Genome evolution and expansion/contraction of gene families

Phylogenetic trees were constructed with 22 single-copy gene families (Supplemental Table S11), which showed that P. ternata and P. pedatisecta were sisters to each other and diverged approximately 11.55 mya. P. pedatisecta was the closest genome to P. ternata (Fig. 2a). Compared with other 13 species indicated that 1,633 gene families were expanded and 4,096 gene families were contracted in P. ternata (Supplemental Table S12). Additionally, 1,112 gene families containing 6,654 genes underwent rapid evolution (Supplemental Table S11). GO and KEGG enrichment analyses were performed and the results revealed that GO analysis enriched in the processes of DNA biosynthesis, RNA-dependent DNA replication and sphingoid biosynthetic and metabolic process (Supplemental Fig. S4a), while KEGG analysis enriched in pathways of Sphingolipid metabolism, Photosynthesis, and Isoquinoline alkaloid biosynthesis (Supplemental Fig. S4b).

Figure 2. 

Evolution and gene family analysis of Pinellia ternata and other representative plant genomes. (a) Phylogenetic analysis of the P. ternata genome and expansions and contractions of gene families. Gene family expansions and contractions are indicated by the numbers in red and blue, respectively. (b) The distribution of shared and specific gene families among six species. Upper diagram showing the distribution of shared gene families among six species. Lower diagram showing gene family complement comparisons among Wolffia australiana, Amorphophallus konjac, Pinellia pedatisecta and Pistia stratiotes, Colocasia esculenta, and P. ternata. There were 8,051 common gene families and 373 P. ternata-specific gene families. (c) Ks distribution among P. ternata and three other species. Different color indicates Ks distribution within and between genomes. (d) Synteny between P. ternata, Amorphophallus konjac and Pinellia pedatisecta. The colored line connects matched gene pairs. The numbers indicate the corresponding chromosomes in each species, Pt, Pp and Ak represents P. ternata, P. pedatisecta and A. konjac, respectively.

The gene family clustering analysis of 14 plant species showed a total of 482,019 genes were clustered into 20,472 gene families (Supplemental Table S11). Among them, a total of 13,200 gene families comprising 28,094 genes were identified in P. ternata genome (Supplemental Table S13). The P. ternata genome assembly was also compared with A. konjac, P. pedatisecta, P. stratiotes, C. esculenta,, and W. australiana. In total, 8,051 gene families common to these six species were identified, and 373 gene families containing 818 genes were specific to P. ternata (Fig. 2b). GO enrichment analysis indicated that the P. ternata-specific genes were mainly enriched in the organic substance biosynthetic process, cellular nitrogen compound biosynthetic process, and cellular biosynthetic process. These enriched genes might be associated with the biosynthetic process of P. ternata (Supplemental Fig. S5a). KEGG enrichment analysis showed that P. ternata-specific genes were mainly enriched for biosynthesis of secondary metabolites, Riboflavin metabolism, Carotenoid biosynthesis, and valine, leucine and isoleucine degradation pathways (Supplemental Fig. S5b).

WGD events and synteny analysis

Whole-genome duplication (WGD) events were investigated in P. ternata, A. konjac, C. esculenta, and P. pedatisecta genomes. P. ternata with the other three species belonging to Aroideae shared one ancient WGD event. The distribution of Ks values in P. ternata showed two peaks at about 0.05 and 0.25 (Fig. 2c), indicating that P. ternata experienced a WGD event twice and that the most recent WGD event occurred independently after P. ternata and P. pedatisecta diverged.

A synteny analysis was conducted with the P. ternata, A. konjac and P. pedatisecta genomes. A total of 13,182 syntenic gene pairs across 379 syntenic blocks between the P. ternata and P. pedatisecta genomes were identified. A total of 11,549 gene pairs across 551 syntenic blocks between the P. ternata and A. konjac genomes were identified. Fragmented conserved synteny was identified for the comparative genome structure between P. ternata and P. pedatisecta (Fig. 2d, Supplemental Fig. S6). In addition, a 1:1 syntenic depth ratio was observed in the P. ternata – P. pedatisecta and P. pedatisecta – P. ternata genomes. For most collinear regions, one chromosome of P. ternata corresponded to one chromosome of P. pedatisecta. For example, Ptchr01, Ptchr02, Ptchr03, Ptchr05, Ptchr06, and Ptchr08 of P. ternata corresponded to Ppchr06, Ppchr03, Ppchr01, Ppchr10, Ppchr04, and Ppchr08 of P. pedatisecta, respectively. The collinearity of P. ternata and P. pedatisecta is obviously better than that of P. ternata and A. konjac (Fig. 2d).

Transcriptome analysis of P. ternata

To gain a better understanding of the medicinal components of P. ternata, a detailed transcriptome analysis was performed on P. ternata tubers of different developmental stages. We compared the tuber of P. ternata in three different developmental stages and that DEGs were identified from O vs Y, and M vs Y comparisons, with 1,394 and 866 DEGs, respectively (Supplemental Table S14). Among the DEGs, 1,039 were specific to O vs Y, and 511 were specific to M vs Y, while 355 DEGs were common in M vs Y and O vs Y (Supplemental Fig. S7).

Next, ClusterGVis was employed to analyze the expression heatmap of the DEGs, and these DEGs among the different developmental stages were clustered into eight expression profiles. A large set of genes (n = 874, clusters 1, 2, 3) expressed at high levels in the Y stage and 386 genes (clusters 1, 2, 3) expressed highly in the M stage, while genes in clusters 6, 7 and 8 (n = 645) expressed at high levels in the O stage (Fig. 3). GO and KEGG analysis were employed to further annotate these genes with high expression levels in different stages. For GO analysis, genes in clusters 1, 2, and 3 were mainly enriched in biological process, including cell cycle process, cell differentiation, single-organism developmental process (Supplemental Fig. S8a). Genes in clusters 4, and 5 had enrichment in biological processes that respond to gibberellin and abiotic stimuli (Supplemental Fig. S8b), while genes in clusters 6, 7, and 8 were mainly enriched in lignin biosynthetic and metabolic processes, phenylpropanoid biosynthetic and metabolic process (Supplemental Fig. S8c). For KEGG analysis, genes with high expression levels in different stages were all enriched in the phenylpropanoid biosynthesis pathways (Fig. 3). Except for enrichment in phenylpropanoid biosynthesis pathways, genes with high expression levels in the Y stage were also enriched in the pathways of pentose and glucuronate interconversions, amino sugar and nucleotide sugar metabolism and starch and sucrose metabolism. Those with high expression levels in the M stage were mainly enriched in tropane, piperidine, and pyridine alkaloid biosynthesis, ascorbate and aldarate metabolism, ubiquinone, and other terpenoid-quinone biosynthesis pathways. The high expression genes in the O stage were mainly enriched in stilbenoid, diarylheptanoid, and gingerol biosynthesis, flavonoid biosynthesis, and plant hormone signal transduction pathways. These combined results revealed that many secondary metabolites are produced during the different growth stages of tubers. Particularly, starch and polysaccharides may be synthesized in the Y stage, and alkaloids are mainly accumulated in the M stage, while stilbenoid, diarylheptanoid and gingerol may be generated with high production in the O stage of tuber.

Figure 3. 

Different gene regulatory landscape in different development stages of Pinellia ternata tuber. (a) Eight expression profile types of DEGs by K-Means clustering, (b) heatmap analysis of all DEGs in eight clusters, (c) pathway enrichment of DEGs in eight clusters.

The KEGG metabolic pathways of DEGs of different clusters were further investigated. Interestingly, genes in the cluster 6 mostly did not express in the Y and M stages but expressed highly in the O stage, and they were enriched in the pathways, including pentose and glucuronate interconversions, plant hormone signal transduction and protein processing in endoplasmic reticulum, suggesting that plant hormones may play a very important role in the tuber growth from the M to O stage.

The TFs in P. ternata were further annotated. A total of 1,564 TFs classified into 68 families were identified, which is approximately 3.5% of the total number of coding genes (Supplemental Table S15). The bHLH family accounted for the largest proportion of coding genes (7.29%), followed by MYB (6.78%), NAC (6.52%), and C2H2 (6.07%). Of all the TFs, 128 were identified as DEGs, which accounts for 8.2% of all TFs and 6.7% of all DEGs. Comparing the three growth stages of tubers, the number of highly expressed TFs were most abundant in the Y stage, followed by O and M stages (Supplemental Fig. S9). The results indicate that these TFs play important roles in regulating tuber growth and development, as well as the bioactive constituent biosysnthesis during this period.

Analysis of key genes involved in ephedrine biosynthesis

The following candidate genes involved in the ephedrine biosynthesis pathway were identified based on the functional annotation and homology search: four phenylalanine ammonia lyase (PAL) genes, two cinnamate: CoA ligase (CNL) genes, two cinnamoyl-CoA hydratase-dehydrogenase (CHD) genes, six 3-hydroxyisobutyryl-CoA hydrolase (CHY) genes, three aldehyde oxidases 4 (AO4) genes, two 3-ketoacyl-CoA thiolase (KAT) genes, three benzaldehyde dehydrogenase (BALDH) genes, five benzoate-CoA ligase (BL) genes, three ThDP-dependent pyruvate decarboxylase (ThPDC) genes, five acetolactate synthase (AHAS) genes, and 13 N-methyltransferase genes (Fig. 4). The candidate genes that were predominantly expressed in tuber compared with the other tissues were as follows: PtPAL PtGene58461, PtCNL PtGene57555, PtCHD PtGene34219, PtKAT PtGene49630, PtBL PtGene38901, PtCHY PtGene12715, PtBALDH PtGene47180, PtAAO4 PtGene00280, and PtAHAS PtGene07361. It is notable that all four PtPAL genes were highly expressed in the root and tuber tissues at different developmental stages. It is possible that the PtPAL genes may be important for allocating the metabolic flux in the underground organs, further contributing to the accumulation of ephedrine in the tuber. For the PtNMTs, phylogenetic analysis, and expression profiles suggested both PtGene15192 and PtGene36030 expressed highly in tuber and had a close genetic relationship with known NMTs from other species (Supplemental Figs S10 & S11). Since NMT is the final key enzyme of ephedrine synthesis, modulating its expression may increase the ephedrine content in P. ternata tubers.

Figure 4. 

Simplified representation of the ephedrine biosynthetic pathway. Top hits for pathway genes identified by blast are highlighted in green. The expression value for each gene is indicated in color on a log2 (FPKM + 1) scale for seven tissue types: leaf, stem, flower, root, young tuber (tuber_Y), medium tuber (tuber_M), and old tuber (tuber_O).

Analysis of candidate genes involved in oxalate biosynthesis and the formation of acrid raphides

The candidate genes involved in oxalate biosynthesis were determined, and the putative oxalate biosynthesis pathways were outlined. Oxalate could be generated either from glyoxylate and oxaloacetate in the tricarboxylic acid (TCA) cycle or from the degradation of ascorbate^[57]. The expression profile revealed the most probable candidate genes: PtACO PtGene58911, PtGLO PtGene39926, and PtMDH PtGene24599 in the TCA cycle and PtAPX PtGene48581 in the ascorbate degradation pathway (Fig. 5).

Figure 5. 

Simplified representation of oxalate biosynthesis and acrid raphides formation. The enzyme names are as follows: glycolate oxidase (GLO); oxaloacetate acetylhydrolase (OXAC); malate synthase (MLS); malate dehydrogenase (MDH); citrate synthase (CTS); aconitase (ACO); isocitrate lyase (ICL); ascorbate peroxidase/ascorbate oxidase (APX/AO); lectin. Top hits for pathway genes identified by blast are highlighted in green. The expression value for each gene is indicated in color on a log2 (FPKM + 1) scale for seven tissue types: leaf, stem, flower, root, tuber_Y, tuber_M, and tuber_O.

Moreover, 14 candidate genes responsible for lectin production were identified, and all were primarily expressed in the tuber (Fig. 5). The expression pattern of lectins may contribute to the formation of acrid raphides in P. ternata. The expression of the candidate genes related to oxalate biosynthesis, together with the lectin genes, may explain the formation of acrid raphides in P. ternata tubers.

Genomically-relevant data is especially limited for studies on the bioactive compounds' biosynthesis mechanism in P. ternata. In the present study, a diploid accession of P. ternata was used for genome sequencing and transcriptome analysis, allowing us to produce a high-quality (LAI > 20) chromosome-level genome assembly of this plant species. Further, the genes involved in ephedrine biosynthesis and acrid raphide formation were identified and found to be highly expressed in the tuber.

As a phenylpropylamino alkaloid, ephedrine is generated in the L-phenylalanine (Phe) pathway, with Phe-derived benzoic acid as an intermediate. The reaction is initiated with a phenylalanine ammonia lyase (PAL) for deamination of phenylalanine to cinnamic acid^[58]. Two possible pathways, i.e., β-oxidative (peroxisome) and non-β-oxidative (cytosol), have been proposed for benzoic acid biosynthesis^[59]. In the β-oxidative routes, four synthases are involved, including the cinnamate: CoA ligase^[60−62], cinnamoyl-CoA hydratase -dehydrogenase^[59,63], 3-ketoacyl-CoA thiolase^[64], and ThDP-dependent pyruvate decarboxylase^[65]. In the non-β-oxidative pathway, two synthases, benzaldehyde dehydrogenase^[66] and acetolactate synthase may be responsible for Phe side-chain shortening for benzoic acid biosynthesis^[67,68]. Moreover, 3-hydroxyisobutyryl-CoA hydrolase is an enzyme that may connect the two pathways by converting trans-cinnamoyl-CoA to benzaldehyde^[69]. Through transcriptome analysis, the genes encoding the above-mentioned biosynthases were obtained in P. ternata^[20]. Nevertheless, none of them were primarily expressed in tubers, failing to explain the accumulation of ephedrine in tubers relative to other organs. In this study, by combining the tissue-specific transcriptome data with genome data, at least one gene for each synthase was found predominantly in the tuber (Fig. 4). Secondary metabolites are usually accumulated and stored in the place of their synthesis^[70]. Thus, the specific expression of the biosynthases is commonly correlated with their accumulation pattern. Unlike Ephedra sinica, which accumulates ephedrine in aboveground parts^[71], P. ternata stores ephedrine in tubers. This indicates that the genes screened in this study are probable as the candidate genes for ephedrine biosynthesis in P. ternata. Moreover, the biosynthesis genes for the conversion of benzoic acid to ephedrine have not been completely discovered in planta, though a transaminase, a reductase, and an N-methyltransferase (NMT) may be involved based on their chemical structure^[17,20] . Through phylogenetic and expression profiles, PtGene15192 and PtGene36030 are considered probable NMTs in P. ternata. Genetically manipulating their expressions may provide insight into breeding P. ternata germplasm with an increased ephedrine content.

In this study, we identified the metabolic pathway of oxalate and further outlined the formation mechanism of acrid raphide in P. ternata. The candidate gene PtGene00795, which was highly expressed in P. ternata tubers, may be important for oxalate accumulation via OXAC (catalyzing oxaloacetate to oxalate)^[72]. There is just one candidate gene each for ICL (PtGene28456) and MLS (PtGene27092), and both have low expression levels in the tuber. This suggests that oxaloacetate rather than glyoxylate is the main precursor of oxalate via the TCA cycle, since oxaloacetate has another source from photosynthesis/glycolysis. In rice, however, glyoxylate is thought to be an efficient precursor for oxalate biosynthesis^[73]. In other plants, both glyoxylate, and oxaloacetate were identified as the precursor^[74]. Ascorbate degradation is another approach for oxalate formation^[75]. Though some catalysis reactions remain clear, it is understood that ascorbate peroxidase/ascorbate oxidase (APX/AO) catalyzes the oxidation of ascorbate^[76]. Among five candidate APX/AO genes, four were highly expressed in tubers (Fig. 5), suggesting that this enzyme is important for the specific biosynthesis of oxalate in tubers. Insecticidal, germicidal, and phlogogenic activities of P. ternata agglutinin (PTA) have been previously reported^[77−79]. As shown in the P. ternata genome data, 14 PTAs were identified and all were highly expressed in tubers, which would facilitate the formation of acrid raphides, as well as the characterization of PTAs.

In conclusion, the genome of P. ternata was sequenced and assembled into 2.08 Gb in 13 pseudo-chromosomes. The expression profile analysis revealed the candidate genes for ephedrine biosynthesis and acrid raphide formation.

Project conception and coordination: Duan Y; experiments design: Xue T; supervision of genome sequencing and assembly: Lin J; genome annotation, synteny and gene family expansion/contraction analyses: Jia H, Wang M; transcriptional data analyses: Zhang Y, Liu X; karyotype analysis: Meng Z; plant materials and DNA extracts maintaining: Chao Q, Zhao F; draft manuscript preparation: Duan Y, Jia H; manuscript revise: Lin J, Xue T, Xue J. All authors reviewed the results and approved the final version of the manuscript

All of the source data, statements of data availability and associated accession codes are available at: https://ngdc.cncb.ac.cn/search/?dbId=&q=PRJCA016404 with Bioproject ID: PRJCA016404. The Pinellia ternata genome assembly and annotation data are available at China National Center for Bioinformation (https://ngdc.cncb.ac.cn/gwh/Assembly/37791/show).