The two-year-old plant of H. japonica in this experiment are in vegetative growth stageused and were collected from the herb garden of Anhui University of Chinese Medicine and identified by Professor Qingshan Yang (Supplementary Figure S1). After washing the fresh H. japonica with ultrapure water, the tissues of leaves (L), roots (R), and stems (S) were separated and then dried with flter paper. All tissues were quickly frozen in liquid nitrogen.

After high-temperature sterilization of the utensils, the leaves, roots, and stems of H. japonica were placed into a mortar and grind them while adding liquid nitrogen. After thorough mixing of the powder, it was placed in a centrifuge tube and the supernatant was collected after centrifugation. The RNA kit (Omega Bio Tek, United States) was used to extract RNA from various tissues. The purity of samples was detected using Thermo NANODROP 2000 ultra micro spectrophotometer, and the concentration and integrity of RNA were detected using Agilent 2,100 bioanalyzer (Supplementary Table S1).

The extracted RNAs were processed using mRNA enrichment and rRNA removal methods. After the obtained mRNAs were fragmented, the first-strand cDNA and second-strand cDNA were synthesized sequentially. The double-stranded cDNA fragments were subjected to end-repair, and then a single ‘A’ nucleotide was added to the 3′ ends of the blunt fragments. The reaction system and program for adaptor ligation were subsequently configured and set up to ligate adaptors with the cDNAs to obtain the cDNA library.

The DNA nanoball sequencing (DNB-seq) platform was used for the sequencing of the leaves, roots, and stems of H. japonica. Single-stranded circle DNA molecules are replicated via rolling cycle amplification, and a DNB which contain multiple copies of DNA is generated. Sufficient quality DNBs are then loaded into patterned nanoarrays using high-intensity DNA nanochip technique and sequenced through combinatorial Probe-Anchor Synthesis. The SOAPNuke (v1.5.2) software (Cock et al., 2010) was used to filter the raw reads obtained from transcriptome sequencing. After removing reads containing adapters, an unknown base (N) content>10%, and of low quality, clean reads were obtained. After clean reads were assembled using Trinity (v2.0.6) software (Haas et al., 2023), transcripts were clustered and deduplicated using CD-HIT (v4.6) software (Fu et al., 2012) to obtain unigenes.

Unigenes were annotated with seven functional databases using hmmscan (v3.0) software (Johnson et al., 2010), Blast (v2.2.23) software (Altschul et al., 1990), and Blast2GO (v2.5.0) software (Conesa et al., 2005): NCBI Nonredundant Protein Sequence Database (NR), NCBI Nucleotide Database (NT), Manually Annotated and Reviewed Protein Sequence Database (SwissProt), Clusters of Eukaryotic Orthologous Groups Database (KOG), Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), Protein Families Database (Pfam). Finally, the functional annotations and classification information corresponding to unigenes were obtained.

Bowtie 2 (v.2.2.5) (Langmead and Salzberg, 2012) and RSEM (v1.2.8) (Li and Dewey, 2011) software were used to calculate the gene expression levels of leaves, roots, and stems of H. japonica, and obtained the standard expression level, expressed as fragments per kill of exon model per million mapped fragments (FPKM). Based on the poisson distribution principle, the differentially expressed genes (DEGs) in leaves, roots, and stems of H. japonica were analyzed and the functions of DEGs were annotated.

The getorf (EMBOSS: 6.5.7.0) and hmmsarch (v3.0) software were used (Rice et al., 2000) to determine the open reading frames (ORFs) of unigenes in H. japonica, and compared the ORFs with the TF protein domains. Unigenes were identified based on the TF family characteristics described in PlantTFDB.

ExPASy (https://web.expasy.org/translate/), MEGA (v5.0) software (Tamura et al., 2011) and CLUSTALX (v1.83) software (Jeanmougin et al., 1998) were used to determine the ORF of SS in H. japonica and compare it with the amino acid sequence of SS in Macleaya cordata, Papaver somniferum, Glycyrrhiza uralensis, Spatholobus suberectus and Glycine max in the NCBI international database. The secondary and tertiary structures of SS in H. japonica were simulated by ESPrip 3.0 (http://esprip.ibcp.fr/ESPri pt/cgi bin/ESPrip.cgi) and Swiss Model (https://swissmodel.epasy.org/), and Pymol (v2.3.2.0) software (Seeliger and de Groot, 2010) were used to visualize its tertiary structure. Phylogenetictrees of SS (CL5764.Contig6) were constructed in MEGA (v5.0) software using the neighbor-joining method with 1,000 bootstrap replicates.

The DNB-seq platform was used for transcriptome sequencing of the leaves, roots, and stems of H. japonica, and the Q30 of high-quality transcriptome reads were not less than 91.66%. After assembly and redundancy removal, 99,404 unigenes were obtained, with total length, average length, N50, N70, N90, and GC content of 158,559,517 bp, 1595 bp, 2335 bp, 1660 bp, 856 bp, and 39.43%, respectively (Supplementary Table S2). Among the obtained unigenes, 59.82% of the unigenes sequences exceeded 1000 bp, and 43.88% of the unigenes sequences exceeded 1500 bp (Supplementary Figure S2). The benchmarking universal single-copy orthologs was used to evaluate the quality of assembled transcripts, 98% of unigenes were successfully matched, indicating good integrity of transcriptome assembly. The RNAsequencing datasets from leaves, roots, and stems of H. japonica were deposited in the NCBI Sequence Read Archive database (accession: PRJNA961922).

Among the 99,404 unigenes, 7,6071 (76.53%), 59,763 (60.12%), 57,396 (57.74%), 59,788 (60.15%), 60,287 (60.65%), 58,349 (58.70%), and 57,404 (57.75%) were annotated by the seven major functional databases of NR, NT, Swissprot, KOG, KEGG, GO, and Pfam, respectively. 31,061 unigenes were annotated by seven major databases, accounting for 31.25%, while 78,989 unigenes were annotated by any one of the seven databases, accounting for 79.46% (Supplementary Table S3). The expression levels of unigenes in 9 samples of leaves, roots, and stems of H. japonica were shown in Figure 1.

An analysis was conducted on 76,071 unigenes annotated in the NR database. The results showed that in the known gene database, the homology relationship between M. cordata (72.28%) and H. japonica was the closest, followed by P. somniferum, Nelumbo nucifera, Aquilegia coerulea and Vitis vinifera (Supplementary Figuer S3). The KOG database annotated 59,788 unigenes, which were classified into 25 categories (Figure 2A). The main annotations were 13,843 unigenes in “General function prediction only”, 47,848 unigenes in “Signal transduction mechanisms”, and 37,403 unigenes in “Function unknown”.

The unigenes annotated in the NR database were further compared to the GO database, resulting in the annotation of 58,349 unigenes, which were divided into three categories: biological process, cellular component, and molecular function (Figure 2B). The biological process were mainly concentrated in the cellular processes (28,823 unigenes) and the metabolic processes (23,941 unigenes). In terms of cellular composition, it is mainly concentrated in the cellular anatomical entity (41,477 unigenes) and the intracellular (18,614 unigenes). The molecular function mainly focus on the binding (39,555 unigenes) and the catalytic activity (38,310 unigenes).

In total, 36,790 DEGs were identified between different tissues of H. japonica. In the comparison of leaves and roots, a total of 31,222 DEGs were identified, of which 14,040 DEGs were up-regulated and 17,182 DEGs were down-regulated in roots. In the comparison of leaves and stems, a total of 14,602 DEGs were identified, of which 6,268 DEGs were up-regulated and 8,334 DEGs were down-regulated in stems. In the comparison of stems and roots, 10,294 DEGs were identified, of which 5,067 DEGs were up-regulated and 5,227 DEGs were down-regulated in stems (Figure 3A). In addition, we identified 1,445 common DEGs among the three comparisons (Figure 3B).

The KEGG database was used for biological functional annotation of DEGs. 11,590 DEGs in the comparison of the leaves and roots were annotated to 134 metabolic pathways (Figure 4A), and 5,589 DEGs in the comparison of the stems and leaves were annotated to 134 metabolic pathways (Figure 4B). These two comparisons were mainly enriched in “Plant-pathogen interaction”, “MAPK signaling pathway-plant”, and “Ribosome”. The 154 DEGs compared between leaves and roots involve the biosynthesis pathway of triterpenoid saponins. 80 DEGs are up-regulated in leaves, while 74 DEGs are down-regulated. The 99 DEGs compared between stems and leaves involve the biosynthesis pathway of triterpenoid saponins. 28 DEGs are up-regulated in the stem, while 71 DEGs are down-regulated. 4,091 DEGs in the comparison of the roots and stems were annotated to 133 metabolic pathways, mainly enriched in “Plant-pathogen interaction”, “MAPK signaling pathway-plant” and “Starch and sucrose metabolism” (Figure 4C). A total of 92 DEGs are involved in the biosynthesis pathway of triterpenoid saponins. 59 DEGs are up-regulated in roots, while 33 DEGs are down-regulated.

60,287 unigenes from the transcriptome of H. japonica were annotated into the KEGG database, including five categories (Figure 5A): cellular processes (2,576 unigenes), environmental information processing (3,693 unigenes), genetic information processing (12,715 unigenes), metabolism (34,168 unigenes), and organismal systems (3,089 unigenes), involving a total of 135 metabolic pathways. In the metabolic of terpenoids and polyketides, unigenes were mainly distributed in “Carotenoid biosynthesis” (Table 1).

Terpenoid backbone biosynthesis (ko00900) and sesquiterpenoid and triterpenoid biosynthesis (ko00909) were two metabolic pathways involved in the biosynthesis of triterpenoid saponins in H. japonica, with a total of 335 unigenes involved. Based on the screening criteria of FPKM>1 of unigenes, there were ten, three, four, three, two, eight, four, and six unigenes for HMGR, MVK, DXS, DXR, MDS, IDI, GPPS, and FPPS in terpenoid skeleton biosynthesis, respectively. In the biosynthetsis pathways of sesquiterpenes and triterpenes, there were two, five, and two unigenes encoding SS, SQLE, and β-AS, respectively (Table 2). The key enzyme unigenes are mostly highly expressed in the leaves and roots of H. japonica, while their expression levels were relatively low in the stems. The relative expression levels of these key enzyme unigenes in leaves, roots and stems were displayed in the form of heat maps (Figure 5B).

Based on the transcriptome database of H. japonica, 2,550 TFs were identified, belonging to 59 TF families (Figure 6A). The most abundant TF families are MYB (300 unigenes), bHLH (210 unigenes), and mTERF (186 unigenes). According to the functional classification of TFs in the KEGG database, 9 TFs were found to be involved in the metabolism of terpenoids and polyketides. Among them, 4 TFs belong to the FHA family, 3 TFs belong to the ABI3VP1 family, 1 TF belongs to the MYB family, and 1 TF belongs to the bHLH family. Among them, CL3296. Contig9, belonging to the bHLH family, was mainly enriched in terpenoid backbone biosynthesis and was closely related to the biosynthesis pathway of triterpenoid saponins in H. japonica. And other TFs were mostly related to the biosynthesis of diterpenes, carotenoids, and secondary metabolites. Protein-protein interaction (PPI) network analysis was performed between these 9 TFs and key enzymes unigenes in the biosynthesis pathway of triterpenoid saponins. The PPI network contained 17 nodes and 107 edges, indicating a high degree of interaction between FPS1, SQE1, and SQE3 with SS (Figure 6B).

In the biosynthesis pathways of sesquiterpenoid and triterpenoid, Squalene is generated by SS catalyzed FPP, and subsequently forms various conformational changes of triterpenoid saponin aglycones. SS is the first key enzyme in this pathway. The ORF length of CL5764. Contig6 of H. japonica was 1278 bp, encoding 425 amino acids. Sequence alignment analysis of SS in H. japonica showed 93.01% sequence homology with SS in M. cordata (McSS, OVA20399.1), 90.48% sequence homology with SS in P. somniferum (PsSS, XP_026393667.1), 88.83% sequence homology with SS in G. uralensis (GuSS, ADG36719.1), 88.83% sequence homology with SS in S. suberectus (SsSS, TKY67615.1), 88.56% sequence homology with SS in G. max (GmSS,NP_001236365.2). The secondary structure of SS in H. japonica was mainly composed of α-helices, with a total of 21 α-helices, as well as 1 TT structure and five η structures (Figure 7A). The three-dimensional structure of SS in H. japonica constructed by Swiss Model has the highest homology with SS from P. somniferum (SMTL ID: A0A4Y7IBE6.1. A), with a similarity of 59%. SS in H. japonica contains six conserved regions (Motif I-VI), with four enzyme activity related sites located within the conserved regions (“VSRSF”, “DTFED”, “DYLED” and “R”) (Figures 7B–D) (Tansey and Shechter, 2001). CL5764. Contig6 encoding SS from H. japonica, as well as SS from different plant species were chosen to construct a phylogenetic tree through multiple sequence alignment (Figure 8).