Multi-Omics Dissection of Huperzine a Biosynthetic Pathway in Huperzia serrata: From Candidate Enzyme Families to Functional Validation of Copper Amine Oxidases

1. Introduction

Huperzia serrata, a traditional Chinese medicinal plant within the family Huperziaceae, is distinguished by its abundance of bioactive compounds, notably the Lycopodium alkaloid huperzine A (HupA) [1,2]. HupA functions as a potent, selective, and reversible inhibitor of acetylcholinesterase and has shown considerable efficacy in the treatment of neurodegenerative disorders such as Alzheimer’s disease [3,4]. Despite its clinical potential, the sustainable supply of HupA presents a significant challenge. The compound is primarily extracted from wild H. serrata plants and a limited number of other Huperziaceae species, which are characterized by extremely slow growth and increasingly scarcity, resulting in a severe supply bottleneck [5]. Therefore, elucidating the complete biosynthetic pathway of HupA is crucial for the development of alternative production platforms through synthetic biology [6].

Previous research has established a foundational understanding of this pathway, tracing its origins to feeding experiments conducted by Ian Spenser and other researcheres [7−9]. These investigations revealed that the biosynthesis of lycopodium alkaloids, including HupA, originates from the decarboxylation of lysine. The initial stages involve intermediates such as cadaverine, ∆¹-piperideine, 4-(2-piperidyl) acetoacetate (4PAA) (or 4-(2piperidyl) acetoacetyl-CoA, 4PAACoA), and pelletierine [7−9]. Further studies identified several critical enzymes, including lysine decarboxylases (LDCs) [10,11], copper amine oxidases (CAOs) [12], and piperidyl ketide synthases (PKSs) [13−15], which are sufficient to convert these primary metabolites. Recently, through D₂O labeling studies and paried transcriptomic and metabolomic analyses, Nett et al. confirmed that de novo HupA biosynthesis is specific to new growth in the shoots of Phlegmariurus tetrastichus and preliminarily elucidated the biosynthesis pathway of HupA, which is broadly divided into four stages [16−18]. As previously menthioned, the pathway initiates with the formation of 8-carbon scaffold precursors, pelletierine [16,18]. In the second stage, pelletierine undergoes sequential modifications, including ketone reduction, O-acetylation, and oxidative desaturation, catalyzed by the short-chain dehydrogenase/reductases (SDRs), the BAHD acyltransferases (ACTs), and the cytochrome P450s (CYP450s), respectively [18]. It is then coupled with 4PAA by a α-carbonic anhydrases (α-CAs) to form the phlegmarane-type scaffold [18]. The subsequent third stage may involve the reconstruction of the core skeleton: the phlegmarane-type scaffold is converted into the lycodane tetracyclic core through a series of steps not yet fully elucidated; this core is then acted upon by isomerases and acetyltransferases, among others, to generate the key intermediate flabellidine [18]. In the final fourth stage, flabellidine undergoes a series of precise oxidative modifications and ring-cleavage reactions mediated by 2-oxoglutarate-dependent dioxygenases (2OGDs) and an α/β-hydrolase (ABH). These reactions convert it successively into HupB, then HupC, and finally yield the target product HupA [16,18] (Figure 1).

Several other studies also have explored the hypothesized biosynthesis pathway of HupA in various Lycopodium species. Comparative analysis of 454-ESTs from H. serrata and P. carinatus have predicted potential genes involved in lycopodium alkaloid biosynthesis, including LDC-likes and some CYP450s [19,20]. Comprehensive transcriptome analysis across all tissues of H. serrata has yielded further insights into the early and late stages of HupA biosynthesis [21]. Recent investigations comparing the in vitro H. serrata thallus with wild H. serrata, utilizing transcriptome and widely targeted metabolome analysis, have identified seven significant enzyme genes and 13 transcription factors (TFs) associated with HupA biosynthesis [6]. Despite these advancements, several aspects of the pathway remain to be elucidated; specifically, conversion process from the phlegmarane scaffold to the lycodane core and the genes corresponding to several hypothesized enzymatic steps require identification and functional validation.

To address these knowledge gaps, we conducted an comparative transcriptomic analysis of various tissues (roots, stems, young leaves, and mature leaves) of H. serrata. Utilizing phylogenetic clustering, gene expression profiling, correlation analysis with HupA content, and protein-protein interaction network prediction, we systematically identified candidate genes potentially involved in HupA biosynthesis. Furthermore, we successfully cloned two candidate CAOs, and demonstrated its crucial roles in amine and alkaloid synthesis and metabolism in planta. Our study not only provides valuable genetic resources but also establishes a foundational framework for fully elucidating the HupA biosynthetic pathway, thereby facilitating its future production through synthetic biology.

2. Materials and Methods

2.1. Plant Materials

The three-branched (THB) H. serrata specimens were collected from Longsheng Huaping Nature Reserve located in Guilin City, Guangxi Zhuang Autonomous Region, China. The roots, stems and leaves of THB H. serrata were either immediately frozen in liquid nitrogen for the purpose of RNA extraction or dried in an oven at 60 °C to assess the contents of HupA and HupB.

2.2. Determination of HupA and HupB Contents

The quantification of HupA and HupB in various tissues of H. serrata was conducted in accordance with the method outlined by Zhu et al. [22], with minal modifications. Specifically, 0.5 g of homogenized powder was precisely weighed and combined with 4 mL of a 0.2% hydrochloric acid solution. The mixture was thoroughly vortexed and subjected to ultrasonic-assisted extraction at 25 °C for 60 min. Following centrifugation, the supernatant was collected, and the residue was re-extracted using the same procedure. The combined supernatants were diluted to a final volume of 8 mL with 0.2% hydrochloric acid, vortex-mixed for homogeneity, and filtered through a 0.22 μm nylon membrane prio to high performance liquid chromatography (HPLC) analysis. HPLC separation was performed on a C18 reversed-phase column (250 × 4.6 mm, 5 μm) maintained at 25 °C. An isocratic mobile phase consisting of acetonitrile and 0.085% formic acid (80:20, v/v) was employed at a flow rate of 1.0 mL/min. Detection was carried out at 309 nm using a VWD detector, with an injection volume was 20 μL.

2.3. RNA Extraction, RNA-seq and Functional Annotation

Total RNA was extracted from various tissues of THB H. serrata, using the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA). RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and RNase free agarose gel electrophoresis. Poly(A) mRNA was enriched using Oligo(dT) beads and fragmented with a fragmentation buffer. These fragments were reverse transcribed into first-strand cDNA using the M-MuLV reverse transcriptase system. Second-strand cDNA was performed using RNase H, DNA polymerase I, dNTP and buffer. The cDNA fragments were then purified with the QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end-repaired, poly(A) were added, and they were ligated to Illumina sequencing adapters. The ligation products were size-selected via agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeqTM 4000 by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China). To ensure high quality clean reads, the original and low-quality data were filtered using fastp (version 0.18.0) [23].

For the full-length transcriptome sequecing, total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA), followed by mRNA enrichment with Oligo(dT) magnetic beads. Enriched mRNA was reverse-transcribed into cDNA using the Clontech SMARTER kit (Takara, Kusatsu, Shiga, Japan). After PCR cycle optimization, double-stranded cDNA was generated; a portion was size-selected (> 4 kb) using BluePippin and mixed with non-size-selected cDNA. SMRTbell libraries were constructed and sequenced on the PacBio Sequel platform. For bioinformatics analysis, raw reads were processed with SMRT Link v5.0.1 [24] to obtain circular consensus sequences, which were classified into full-length non-chimeric (FLNC) reads. FLNC reads were clustered by Iterative Clustering for Error Correction software, polished with Quiver (accuracy ≥ 99%), and low-quality isoforms were error-corrected using Illumina short reads. Redundant sequences were removed with software CD-HIT-v4.6.7 (0.99 identity). Clean reads from the second-generation transcriptome were mapped to the full-length transcriptome and quantified using RSEM software [25].

All assembled sequences were subjected to BLAST analysis across various databases, including kyoto encyclopedia of genes and genomes (KEGG), the NCBI nonredundant protein sequences (Nr), a manually annotated and reviewed protein sequence database (SwissPro), gene ontology (GO), cluster of orthologous groups of proteins (COG), eukaryotic ortholog groups (KOG), and the translation of EMBL (TrEMB) [26]. The GO functional annotation was derived from the annotations of transcripts in Swissprot and TrEMBL. The KEGG Automatic Annotation Server was employed to obtain KEGG annotation information. Open reading frmes were predicted by ANGEL [27], protein domains by Pfam [28] and SMART, transcription factors by PlantTFdb ((http://planttfdb.cbi.pku.edu.cn/), R-genes by PRGdb (http://prgdb.crg.eu/wiki/Main_Page), and alternative splicing by Cogent [29].

2.4. Differential Gene Analysis

The abundance of each gene was normalized to Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) score. Genes with a false discovery rate (FDR) below 0.05 and an absolute fold change greater than one were identified as differentially expressed genes (DEGs) [30]. Principal Component Analysis (PCA) was employed to elucidate the relationships among the samples and was conducted using R package models (http://www.r-project.org/) in this study. Furthermore, a functional enrichment analysis using the KEGG database was performed to determine which DEGs were significantly enriched in KEGG metabolic pathways.

2.5. Identification of Candidate Genes Involved in the HupA Biosynthetic Pathway

This study employed two distinct methodologies for the systematic identification and selection of candidate genes potentially involved in HupA biosynthesis.

2.5.1. Candidate Gene Screening Combining Phylogenetic and Co-Expression Analysis

The screening process commenced with the retrieval of Hidden Markov Models (HMMs) corresponding to nine protein families previously identified as participatants in the HupA biosynthetic pathway in P. tetrastichus (Table 1) [16,18]. These HMM profiles were sourced from the Arabidopsis database (TAIR, https://www.arabidopsis.org/) and utilized for domain-based screening of the H. serrata transcriptome-encoded proteins (Figure S1). The resultant sequences underwent multiple sequence alignment with functionally validated homologous proteins, followed by the construction of a phylogenetic tree to preliminarily assess potential functions based on evolutionary clustering patterns (Figure S1). Subsequently, by integrating tissue-specific expression profiles, genes exhibiting high expression levels in young leaves were selected. The Pearson correlation coefficients between their expression levels and the contents of HupA and HupB were calculated. Genes demonstrating statistically significant positive correlations (p < 0.05) were retained as high-confidence candidates (Figure S1).

2.5.2. Candidate Gene Screening Leveraging Differential Expression and Functional Annotation

The alternative approach commenced with the identification of transcripts significantly upregulated in young leaves compared to mature leaves across all genes expressed in H. serrata (Figure S1). Subsequently, a correlation analysis was conducted between their expression levels and HupA content. Genes exhibiting strong positive correlations (correlation coefficient > 0.8) were selected for further analysis (Figure S1). Following functional domain annotation, these genes were mapped to protein families with experimentally validated roles in prior studies, facilitating the inference of their potential catalytic or regulatory functions within the HupA biosynthetic pathway (Figure S1).

2.5.3. Phylogenetic and Correlation Analysis

Multiple amino acid sequence alignments were conducted using MAFFT (v7) [40]. Phylogenetic trees were generated employing the Neighbor-Joining method in MEGA (v11) [41], with 1000 bootstrap replications. Pearson correlation analysis was performed to investigate the relationships between the expression of alkaloid-related genes and the relative content of HupA and HupB. Multiple group comparisons were analyzed using one-way analysis of variance, followed by Tukey’s multiple comparisons test, implemented in the OmicShare online tool (http://www.omicshare.com/tools) [42]. A p-value less than 0.05 was considered statistically significant.

2.6. Network Analysis of Key Alkaloid Biosynthesis Genes

The significant correlation network, characterized by a p-value of less than 0.01 and a correlation coefficient exceeding 0.9, or an interaction score greater than 700, was visualized using Cytoscape software [43].

2.7. Gene Cloning and Agrobacterium-Mediated Transient Expression

Total RNA was extracted from a composite sample comprising roots, stems, young leaves, and mature leaves of H. serrata using the Plant RNA Prep Pure Plant Kit (Vazyme, Nanjing, China) in accordance with the manufacturer’s protocols. The quantity, purity, and integrity of the RNA were evaluated through 1% agarose gel electrophoresis and spectrophotometry (Thermo Fisher Scientific, Waltham, Massachusetts, USA). cDNA synthesis was conducted utilizing the TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Vazyme, Nanjing, China). The coding sequence (CDS) of HsCAO1 was amplified using 2×SparkHiFi Max Master Mix (SparkJade, Jinan, China) along with gene specfic primers (Table S1). The cDNA of HsCAO genes was determined based on the CDS of HsCAO1 employing a SMARTer RACE 5’/3’ kit (Takara, Kusatsu, Shiga, Japan), with specific primers (Table S18). The CDSs of HsCAO1, HsCAO2, and HsCAO3 were lingated into the overexpressing vector, respectively. The empty vector (EV) or recombinant vector was introduced into the leaves of N. benthamiana via Agrobacterium-mediated transformation. Gene expression levels were assessed using DNA or cDNA templates from transgenic and EV tobacco leaves 48 to 72 h post-injection.

2.8. Determination of Amines and Alkaloids in the Leaves of N. benthamiana

The quantification of the three biogenic amines—putrescine, spermidine, and spermine—was conducted using reversed-phase HPLC with pre-column derivatization. The samples underwent extraction with perchloric acid, followed by ultrasonication and centrifugation. The resulting supernatant was derivatized through the addition of sodium hydroxide, saturated sodium bicarbonate, and dansyl chloride in acetone, and subsequently incubated at 50 °C in the dark for 40 min. Ammonia was then introduced, and the mixture was allowed to stand at room temperature for 30 min prior to dilution with methanol and filtration. Separation was achieved using a Compass C18(2) column (250 mm × 4.6 mm, 5 μm) maintained at 30 °C, with a mobile phase gradient of methanol and water at a flow rate of 1 mL/min, and detection at 254 nm. For quantification purposes, standard solutions of each amine were derivatized in an identical manner and injected under the same conditions; calibration curves were generated by plotting peak area against content.

The quantification of nicotine was conducted using reversed-phase HPLC. The sample was initially ground, and 0.2 g was extracted with 10 mL of 75% methanol containing 0.1 mol/L NaOH through ultrasonication for 1 h, followed by centrifugation and filtration. The separation process was executed on a Compass C18(2) column (250 mm × 4.6 mm, 5 μm) maintained at 30 °C, utilizing a mobile phase composed of methanol and 0.1% triethylamine (45:55) at a flow rate of 1.0 mL/min, with detection occurring at 260 nm. Quantification was accomplished via a standard curve generated using nicotine standards dissolved in methanol.

The quantification of total alkaloids was conducted via a colorimetric method employing bromocresol green. Dried and powdered samples (0.1 g, 80-mesh) were subjected to extraction with 1 mL of 80% ethanol through ultrasonication for 60 min, followed by centrifugation at 8000 g for 10 min. The resulting supernatant was treated with reagent 1 and reagent 2, allowed to react for 5 min, and subsequently subjected to chloroform extraction. After a shaking period of 40 min, the chloroform layer was analyzed at a wavelength of 416 nm. Quantification was achieved using a standard curve defined by the equation y = 2.7614x + 0.0116.

3. Results

3.1. HupA and HupB Content in Different Tissues of H. serrata

H. Serrata (Figure 2A-B, Figure S2A-C), a terrestrial species thriving in shaded, damp, mossy, and acid soil beneath woodlands, exhibits a notably slow growth rate, typically requiring 15 to 20 years to reach maturity from spore germination [2]. As the plant ages, it develops taller and more dichotomous branches (Figure S2B-D). Previous studies utilizing P. tetrastichus have demonstrated that HupA accumulates to a significant level across all tissues, with de novo HupA biosynthesis being specific to new growth in the shoots of P. tetrastichus [16,17]. In this study, we quantified the contents of HupA and HupB in the roots, stems, young leaves, and mature leaves of H. serrata with varying numbers of branches using HPLC. The results indicated that HupA and HupB were predominantly accumulated in the young leaves of H. Serrata, irrespective of the number of branches (Figure 2C-D, Figure S2E-F). This finding aligns with the accumulation pattern of HupA observed in the various tissues of P. tetrastichus [16,17], suggesting that active biosynthesis of HupA also occurs in the newly developed tissues of H. serrata.

3.2. Analysis of Transcriptomic Differences Between Various Tissues of H. serrata

To further elucidate the molecular mechanisms underlying HupA biosynthesis in H. serrata, transcriptome sequencing was conducted on variable tissues, including roots, stems, mature leaves and young leaves of THB H. serrata, utilizing Illumina HiSeqTM 4000 platform. The transcriptome sequencing of these tissues of H. serrata yielded between 5.92 and 7.36 Gb of raw bases, with Q20 and Q30 values exceeding 98.19% and 94.68%, respectively, and a GC content ranging from 43.48% to 44.05% (Table S2). Following quality control analysis, 5.87 to 7.29 Gb of clean bases were obtained, comprising 39,502,296 to 49,043,080 reads (Table S2). The average GC content of the clean bases was determined to be 43.39% and 43.97%, with Q20 and Q30 above 98.33% and 94.85%, respectively (Table S2), indicating high-quality data suitable for subsequent bioinformatic analysis.

In total, the analysis of these transcription profiles targeted the expression intensity of 51,804 genes. PCA of global gene expression effectively distinguished the variances among different samples, while samples within the same group clustered closely, indicating strong biological repeatability within the sample group and high reliability of the transcriptome data (Figure 3A,Table S3).

To identify DEGs across various tissues, we calculated the log2-fold changes in gene expression levels among the six tissues. Genes exhibiting an absolute fold change greater than 1 and a FDR less than 0.05 following correction were classified as DEGs. A total of 17,259 DEGs were identified between young leaves and other tissues of the THB H. serrata (Table S4). Furthermore, 1,953 DEGs were common to all compared groups (Figure 3B, Table S5).

KEGG pathway classification was conducted to categorize the functions of 1,953 DEGs in the young leaves of the THB H. serrata (Table S6). A total of 417 candidate genes with pathway annotations were assigned to 81 KEGG pathways (Table S6). Among these, 156 genes were classified under the “Biosynthesis of other secondary metabolites (ko01110)” pathway, which includes genes such as chalcone synthase, shikimate O-hydroxycinnamoyltransferase, and ornithine decarboxylase (Table S7).

To investigate comprehensive and accurate transcript isoforms, including alternative splicing, novel genes, and fusion transcripts, and to significantly enhance the precision of quantification and isoform-level analysis in short-read RNA-seq data from individual tissues, full-length transcriptome sequencing of a mixed sample comprising six tissues was conducted on the PacBio Sequel II platform. A total of 35.63 Gb of raw reads were generated, encompassing 20,659,204 subreads with an average length of 1,725 bp (Table S8). Following reads clustering and correction, 41,414 polished high-quality isoforms were obtained. To further refine the accuracy of PacBio ISO-seq results, Illumina RNA-seq results were employed for isoforom correction. Consequently, a total of 34,514 non-redundant transcripts with an N50 of 2,257 bp, a length distribution ranging from 61 to 8,058 bp and GC content of 43.74% (Figure 4A, Table S9).

Within the entire set of transcripts,, 32,858 (95.2%), 31,607 (91.6%), 26,845 (77.8%) and 20,646 (59.8%) unigenes were successfully annotated in NR, KEGG, SwissProt and KOG databases, respectively (Figure 4B). In total, 32,900 unigenes were annotated in at least one database, representing approximately 95.32% of all genes (Figure 4B). In NR database, 95.2% of the unigenes exhibited similarity with other plant species, including 31.1% (10,731 unigenes) of all ungenes with Selaginella moellendorffii, 18.3% (6,330 unigenes) with Physcomitrella patens, 1.8% (618 unigenes) with Cinnamomum micranthum f. kanehirae, 1.5% (510 unigenes) with Amborella trichopoda, and 1.5% (504 unigenes) with Quercus suber (Figure 4C, Table S10). Furthermore, 13,557 genes were categorized into 25 groups in the KOG database, including “general function prediction only” (4,471 genes), “posttranslational modification, protein turnover, chaperones” (3,304 genes), and “signal transduction mechanisms” (2,501 genes) (Figure 4D, Table S11). TFs are pivotal in plant development and metabolite biosynthesis processes by modulating gene expression. In this study, 1,201 putative TFs, classified into 39 families, were identified (Figure 4E, Table S12). The majority of TFs belonged to the ERF family (160 genes, 13.32%), followed by the WRKY (99 genes, 8.24%), GRAS (94 genes, 7.83%), Trihelix (91 gnes, 7.58%), and bZIP (88 genes, 7.33%) (Figure 4E, Table S12). Additionally, a total of 749 alternative splicing events were predicted. The analysis revealed that intron retention was the most frequent event, occurring 372 times (49.67%), followed by alternative 5’ donor sites occurring 158 times (21.11%), alternative 3’ donor sites occurring 154 times (20.56%), while alternative last exons were the least frequent, with only 10 occurrences (1.34%) (Figure 4F).

3.3. Candidate Genes Involved in the Biosynthetic Pathway of HupA

Through the application of phylogenetic tree clustering, expression level analysis, and correlation analysis between gene expression and metabolite content, we identified 71 genes potentially involved in the backbone synthesis of HupA. These genes belongs to multiple enzyme families, including LDC, CAO, chalcone synthase (CSH), neomenthol dehydrogenase-like (NMRL), CYP450, BAHD ACT, and α-CA (Figure 5, Figure S3-S9). Furthermore, 28 genes likely participating in the modification of the HupA skeleton were identified, primarily encompassing 2OGD and caffeoyl shikimate esterase (CSE) (Figure 6, Figure S9-S10). Among these, five CSH genes (HsCSH-2, HsCSH-3, HsCSH-4, HsCSH-10, HsCSH-11), five CYP450 genes (HsF3’M-1 to HsF3’M-5), five BAHD ACT genes (HsACT-1, HsACT-2, HsHCT-1, HsHCT-2, HsHCT-3), thirteen α-CA genes, four 2OGD genes (Hs2OGD-1, Hs2OGD-2, Hs2OGD-3, HsGA2OX-3), and one CSE gene (HsCSE-15) were highly expressed in young leaves and demonstrated a significant correlation with HupA content.

To mitigate the confounding effects of tissue-specific gene expression on the analysis of candidate differential genes involved in HupA biosynthesis, a differential gene expression analysis conducted between young and mature leaves, with respect to HupA content. This analysis identified 3,801 up-regulated genes in young leaves (Figure 7A-C, Table S13). Among these, 113 genes exhibited a significant correlation with HupA content (Pearson correlation coefficient (PCC) > 0.8), encompassing 422 enzyme families (Table S14). Notably, this included 86 genes from seven protein families implicated in HupA biosynthesis: LDC, CSE, NMR, CYP450, α-CA, 2OGD, and CSE, all of which demonstrated high expression levels in young leaves (Figure 7D, Table S15).

3.4. Network Analysis of Key Alkaloid Biosynthesis Genes

To further investigate additional genes implicated in the regulation of HupA biosynthesis, the previously identified seven candidate genes were analyzed using protein-protein interaction (PPI) network analysis. This analysis demonstrated that these proteins not only interact with each other but also associate with multiple proteins involved in HupA biosynthesis. These include shikimate O-hydroxycinnamoyltransferase and acyltransferase from the BAHD acyltransferase family, gibberellin 2-beta-dioxygenase from the Hs2OGD family, and the CYP782B7 subfamily protein from the CYP450 family (Figure 8A, Table S16). Furthermore, interactions were identified with calcium-transporting ATPase, cyanate hydratase, caffeoyl-CoA O-methyltransferase, cinnamoyl-CoA reductase, trans-cinnamate 4-monooxygenase, 4-coumarate-CoA ligase, flavonoid 3’,5’-hydroxylase, and zeaxanthin epoxidase, among others (Figure 8B, Table S17).

Furthermore, the seven candidate genes previously identified for HupA biosynthesis were utilized for TF interaction network analysis. The findings revealed that the proteins encoded by these genes interact with members of various TF families, including MYB, NF-YC, GRAS, ERF, bHLH, GRF, and SAP (Figure 8C, Table S18).

3.5. HsCAOs Was Involved in the Alkaloid Synthesis and Metabolism in Nicotiana benthamiana

As previously mentioned, CAOs have the potential to catalyze cadaverine to form ∆¹-piperideine, a crucial precursor in the initial stage of HupA biosynthesis [1,5,16,18]. Sun et al. successfully cloned a CAO (HsCAO) from the roots of H. serrata [12], reporting that HsCAO could utilize not only cadaverine as a substrate but also several others, including putrescine, spermidine, histamine, benzylamine, tyramine, 4-aminomethylpyridine, 3- aminomethylpyridine, and 2-aminoethylpyridine [12]. In this study, we identified no isoforms identical to HsCAO (herein referred to as HsCAO1) in the transcriptome and full-transcriptome data. Consequently, primers were designed based on the CDS of HsCAO1, facilitating the successful cloning of two additional highly conserved sequences through the 3’- and 5’-RACE methods,, which we have designated as HsCAO2 and HsCAO3, respectively. To evaluate the putative substrates of these CAOs, we attempted to express and purify the three recombinant HsCAO proteins using various prokaryotic protein expression systems and the Pichia pastoris expression system. Regrettably, no active HsCAO proteins were obtained, including HsCAO1. Subsequently, we overexpressed HsCAOs and assessed the content of several amines and alkaloids in the leaves of N. benthamiana. The results indicated a decrease in putrescine content and a significant increase in nicotine content in HsCAO-overexpressed lines (HsCAO-OX) compared to those in the empty vector-injected N. benthamiana leaves. Furthermore, the content of spermidine, spermine, and total alkaloids in HsCAO-OX varied compared to those in EV, suggesting the complex roles of different HsCAOs in the amines and alkaloid synthesis and metabolism in planta.

Figure 9. The contents of amines and alkaloids in HsCAO-overexpressed N. benthamiana lines. (A) The schematic representation of the overexpression vector structure. (B) The schematic representation of the tobacco injection procedure. (C)~(I) The contents of putrescine, spermine, spermidine, nicotine, and total alkaloids in HsCAO-overexpressed N. benthamiana lines, respectively. M: DL2000 maker; WT: the wild type N. benthamiana; V: N. benthamiana injected with EV; 1,2, and 3: N. benthamiana injected with HsCAO1-OX, HsCAO2-OX, and HsCAO3-OX vectors, respectively.

Figure 9. The contents of amines and alkaloids in HsCAO-overexpressed N. benthamiana lines. (A) The schematic representation of the overexpression vector structure. (B) The schematic representation of the tobacco injection procedure. (C)~(I) The contents of putrescine, spermine, spermidine, nicotine, and total alkaloids in HsCAO-overexpressed N. benthamiana lines, respectively. M: DL2000 maker; WT: the wild type N. benthamiana; V: N. benthamiana injected with EV; 1,2, and 3: N. benthamiana injected with HsCAO1-OX, HsCAO2-OX, and HsCAO3-OX vectors, respectively.

4. Discussion

The biosynthesis of natural products via biotechnological systems constitutes a progressive approach [44]. In this context, synthetic biology has facilitated new pathways for the efficient large-scale production of bioactive compounds [44]. As an emerging and rapidly advancing interdisciplinary field within biotechnology, synthetic biology has been predominantly utilized with microbes and plants for the synthesis of bioactive compounds [45]. This methodology is particularly advantageous for bioactive compounds with limited production, such as HupA, due to the scarcity of plant resources, the minimal quantities present in plants, and the complexity of their chemical synthesis [46].

The identification of key missing enzymes, and consequently the elucidation of the biosynthetic pathway of bioactive compounds in medicinal plants through proteomics, metabolics, particularly transcriptomics and genomics, forms the basis for the heterologous biosynthesis of these bioactive compounds [46]. Unfortunately, extant lycophytes, especially homosporous Lycopodiaceae, inculding H. serrata and and other HupA-producing species, possess large genome size, high chromosome numbers, and significant chromosome polyploidization [47,48], which impede the acquisition of high-quality genome information for these species. To date, only a limited number of nuclear and organellar genomes of Lycopodiaceae species have been characterized [47,49,50]. As an alternative, identification of key missing enzymes via transcriptomics may present a viable strategy. This approach has been successfully employed in the elucidation of biosynthetic pathway for various natural products, including colchicine [51], strychnine [52], galanthamine [53], thereby confirming the feasibility of analyzing biosynthetic pathways based on transcriptome analysis.

By integrating PacBio long-read and Illumina short-read sequencing technologies, this study successfully assembled a high-quality transcriptome database for H. serrata. This dataset not only serves as a valuable resource for research on HupA biosynthesis but also offers new insights into the evolutionary biology of the Lycopodiaceae family. Comparisons with the NR database revealed that H. serrata exhibits the highest sequence similarity to S. moellendorffii and P. patens (Figure 4C), further supporting its primitive status in the evolutionary tree of land plants. KOG functional classification indicated that “General function prediction only,” “Posttranslational modification, protein turnover, chaperones,” and “Signal transduction mechanisms” were the categories with the highest number of genes (Figure 4D). This pattern reflects the biological characteristics of H. serrata as a perennial slow-growing plant [54]. Its extended life cycle and slow growth rate may necessitate sophisticated mechanisms for maintaining protein homeostasis to cope with long-term environmental stresses. Meanwhile, a complex signal transduction network likely enables it to precisely perceive and respond to specific environmental conditions such as shaded habitats and acidic soils [55].

The study determined that the contents of HupA and its precursor HupB were markedly elevated in young leaves compared to other tissues of H. serrata (Figure 2C-D, Figure S2E-F). This finding suggests that young leaves serve as the primary site for HupA biosynthesis, which was consistent with the finding observed in P. tetrastichus [17], thereby offering a definitive target for gene screening within this pathway. Through phylogenetic clustering, expression level analysis, and correlation analysis between gene expression and metabolite content, 71 candidate genes potentially involved in the skeleton synthesis of HupA and 28 candidate genes potentially involved in its modification were identified (Figure 5A-G, Figure 6A-B). An association analysis between DEGs in young and mature leaves and the contents of HupA and HupB identified 113 genes significantly correlated with HupA content (PCC > 0.8) (Table S13). The partial overlap between the gene sets obtained through the two screening methods underscores the reliability of the strategy and provides essential entry points for elucidating this complex pathway.

In the initial stage of HupA biosynthetic pathway, three catagories of enzymes—LDC, CAO, and PKS—are involved (Figure 1). LDC, which catalyzes the decarboxylation of L-lysine to yield cadaverine, is identified as the initial and rate-limiting enzyme in the HupA biosynthetic pathway [6,10]. To date, six HsLDCs have been successfully cloned and reported; however, only two possess the capacity to decarboxylate L-lysine to yield cadaverine [10,11,56]. Additionally, these two HsLDCs could decarboxylate L-ornithine to produce putrescine [10,56]. Phylogenetic clustering analysis in our study indicated that these two HsLDCs are classified into one group (Figure 5A). Furthermore, 51 isoforms potentially encoding LDCs were identified in our study (Figure S3), suggesting that these LDCs may be generally involved in the synthesis of various alkaloids and may have developed substrate specificity and subcellular localization specificity through gene duplication and functional differentiation [57]. Comparative transcriptome analysis has identified nine HsLDC isozyme genes highly expressed in young leaves, among which HsLDC-9 exhibited the highest expression level in young leaves, with its expression pattern showing a strong positive correlation with HupA content (Figure 5A). It is hypothesized that HsLDC-9 might be localized in the cytoplasm and responsible for providing the cadaverine precursor for HupA synthesis. Further investigation should be conducted to verify this hypothesis.

CAO is the subsequent key enzyme downstream of LDC in the HupA biosynthetic pathway, catalyzing the oxidative deamination of cadaverine to produce Δ¹-piperideine [12]. This reaction introduces a critical aldehyde group and constructs the core piperidine ring structure of HupA [58]. Seven isoforms potentially encoding CAOs were identified in our study (Figure S3); however, the expression of all these isoforms in young leaves does not strongly correlate with HupA content (Figure 6B). As HsCAO1 was reported to be a functional CAO involved in HupA biosynthesis [12], we successfully cloned it. Interestingly, other two highly conserved CDSs were also identified. However, we were unable to purify the recombinant HsCAOs because of their insoluble nature. Alternatively, we overexpressed them in the leaves of N. benthamiana, respectively, and the results demonstrated their different roles in amines and alkaloids synthesis and metabolites (Figrue 9E-I). Bewilderingly, Nett et al. reported the HsCAO1 could not catalyze cadaverine to be Δ¹-piperideine in vitro [16], which cotrodicts the result reported by Sun et al. [12]. Phylogenetic clustering analysis in our study indicated that HsCAOs and functional CAOs identified from P. tetrastichus were not closely related (Figure 6B), indicating their functional differentiation. Taken together, the function of HsCAOs requires further investigation.

The CYP450 enzyme superfamily represents one of the most catalytically diverse classes of enzymes in plants, capable of facilitating reactions such as hydroxylation, epoxidation, dealkylation, and desaturation [59]. Prior research on P. tetrastichus has revealed that the CYP450 enzyme PtCYP782C1 is involved in the desaturation reaction that forms the HupA skeleton [18]. The CYP782B7 subfamily members identified in this study are linked to the synthesis of specialized alkaloids in other plants, and their elevated expression in young leaves aligns with their established roles in alkaloid diversification (Figure 5E). This indicates that these enzymes may play a crucial role in HupA biosynthesis. The BAHD ACT is noted for its extensive substrate specificity and catalytic versatility [60]. In P. tetrastichus, PtACT-1 has been verified to catalyze the formation of the amide bond within the HupA skeleton [18]. This investigation identified five BAHD family enzymes (e.g., HsHCT-1, HsACT-2) (Figure 6F), which may exhibit analogous functions.

CAs are typically involved in CO₂ fixation and pH balance regulation [61]. However, recent research has demonstrated that certain members of this family can directly participate in the biosynthesis of lycopodium alkaloids [18]. For example, PtCAL-3 catalyzes the stereoselective condensation of Δ¹-piperideine and 3-oxoglutaric acid to produce (S)-4-(2-piperidyl)acetoacetic acid, while PtCAL-1a/1b and PtCAL-2a/2b must be co-expressed to catalyze subsequent condensation reactions [18]. This study identified multiple carbonic anhydrase genes highly homologous to PtCAL (Figure 5G), with expression levels strongly positively correlated with HupA and HupB content, suggesting their potential involvement in similar condensation reactions.

2OGDs are frequently involved in the hydroxylation of secondary metabolites in plants [62]. In P. tetrastichus, five 2OGD enzymes have been shown to participate in the modification of the HupA skeleton, including pyridone ring formation, piperidine ring cleavage, and A-ring carbonylation [16,18]. In this study, Hs2OGD-1, Hs2OGD-2, and Hs2OGD-3 were highly homologous to Pt2OGD-1 and significantly positively correlated with HupA and HupB content (Figure 6A). It is hypothesized that these enzymes may catalyze the cleavage of the piperidine ring of HupB, resulting in the loss of one carbon atom and the generation of HupC.

The biosynthesis of HupA is likely governed by a complex, multi-tiered network that integrates developmental signals, environmental factors, and transcriptional regulation. PPI analysis has indicated that the candidate proteins not only interact with each other but also exhibit extensive connections with enzymes involved in phenylpropanoid metabolism (e.g., C4H, 4CL, CCR) (Figure 8A-B). This suggests potential crosstalk between HupA biosynthesis and pathways such as phenylpropanoid metabolism and plant hormone signal transduction. Such metabolic network interconnections may represent a common strategy for plants to coordinate the synthesis of multiple defense compounds. TF interaction analysis revealed that multiple candidate genes could interact with members of the MYB, bHLH, and ERF TF families (Figure 8C). These factors play central roles in the regulation of plant secondary metabolism: MYB transcription factors often act as hubs integrating developmental and environmental signals [63]; bHLH transcription factors frequently form complexes with MYB to activate metabolic pathways; ERF transcription factors typically respond to hormonal signals such as jasmonic acid and ethylene to coordinate defense responses [64,65]. For example, the MYB-bHLH-WD40 complex regulates flavonoids and anthocyanin synthesis in other plants [66,67], and similar mechanisms may exist in the regulatory network of HupA synthesis.

In summary, this study presents the inaugural comprehensive multi-tissue transcriptomic atlas of H. serrata, serving as a valuable resource for molecular research on this species. Utilizing a multi-dimensional screening strategy that integrates expression patterns, co-expression analysis, and evolutionary relationships, the study successfully identified multiple high-confidence candidate genes and constructed a preliminary regulatory network. Nonetheless, this study has certain limitations: A majority of functional predictions are based on bioinformatics analysis and lack validation through in vivo or in vitro experiments; the complete biosynthetic pathway of HupA, particularly the isomerization steps of key intermediates, remains to be elucidated; and the current data do not fully explain the low accumulation of HupA in H. serrata or its ecological function and evolutionary origin. Future research should focus on systematically validating the roles of candidate genes through functional experiments to further refine the biosynthetic pathway and regulatory mechanisms of this important medicinal compound.

5. Conclusions

This study employed integrated multi-omics analysis to systematically elucidate the tissue distribution pattern of HupA in H. serrata. It constructed a high-quality transcriptomic resource, identified multiple high-confidence candidate genes, and preliminarily uncovered their potential regulatory network. These findings not only establish a robust foundation for fully elucidating the biosynthetic pathway of HupA but also provide critical targets for enhancing its production through synthetic biology strategies. Future research should focus on the functional validation of the candidate genes and an in-depth examination of the regulatory mechanisms, ultimately achieving sustainable production of HupA to meet clinical demands while conserving this valuable medicinal plant resource.