Individual B. balsamifera plants were selected for vegetative propagation at the germplasm repository of the tropical medicinal plant greenhouse (19°30′N, 109°20′E; temperature, 28°C; humidity, 60%; light application time, 9 h) of Ministry of Agriculture in Danzhou city, Hainan province, which is linked with the Tropical Crops Genetic Resources Institute, Chinese Academy of Tropical Agricultural Sciences. The nursery garden has red clay soil, partly sandy. Every plant was grown for 3 months. The identical plants were separated into five groups based on leaf count and phenotype, each with 20 plants. Each MeJA solution was produced in 1 L of ddH₂O, which provided sufficient volume for foliar applications. Each plant was continuously sprayed with 20 mL of a surface-treating chemical. In addition to the control (distal water), each foliar treatment received equal exposure to MeJA doses of 1 mmol/L and 10 mmol/L. A volume of 20 ml of surface-treating chemicals were added to the spray bottle and equably sprayed at 20 cm from the plants. Each treatment was given three times in total. After the treatment, mature leaves were gathered at 24 h, 48 h, and 120 h. The samples were immediately stored in liquid nitrogen and kept at −80°C, until they were used for qPCR, RNA sequencing (RNA-seq), and metabolite analysis.

The collected B. balsamifera leaves were immediately placed in a mortar and ground with liquid nitrogen, and ground into a powder. We accurately weighed 2 g of the ground powder into a 50-mL centrifuge tube and added 25 mL of ethyl acetate, operating at 40 kHz and 400 W. We performed an ultrasonic extraction for 30 min, left it standing, and filtered it. Add 1 mL of the filtrate to a 10-mL volumetric flask, along with 1 mL of methyl salicylate internal standard solution. Next, fix the volume with ethyl acetate and filter through a 0.22-μm microporous filter membrane. The content of L-borneol was determined by GC-MS, HP-5 quartz capillary column (0.32 mm × 30 m, 0.25 μm), starting at 80°C for 2 min, then heating at 5°C/min to 100°C, 20°C/min to 200°C. The temperatures of the inlet and FD detector were 220°C and 240°C, respectively. The injection volume is 0.6 μL, without shunt.

TRIzol reagent (Invitrogen) was used to extract total RNA, and the manufacturer’s instructions were followed (El-Sappah et al., 2023). Using the NanoPhotometer1 spectrophotometer (IMPLEN, CA, USA) and the Qubit1 RNA Assay Kit in Qubit1 2.0 Flurometer (Life Technologies, CA, USA), total RNA purity and concentration were assessed. RNA integrity was evaluated using the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) and the RNA Nano 6000 Assay Kit. The RNA sample preparations employed a total of 1.5 g of RNA as input material for each sample. Utilizing the NEBNext1 UltraTM RNA Library Prep Kit for Illumina1 (NEB, USA) in accordance with the manufacturer’s instructions, sequencing libraries were created, and index codes were applied to assign sequences to specific samples. The TruSeq PE Cluster Kit v3-cBot-HS (Illumina) was used on the acBot Cluster Generation System to cluster the index coded sample data in accordance with the manufacturer’s recommendations. Using an Illumina Hiseq 4000 platform, the library preparations were sequenced after cluster creation, and paired-end reads were produced. Fastq format raw readings were processed using internal perl scripts as part of the quality control stage. By eliminating adapter- and ploy-N-containing reads and low-quality reads from the raw data in this stage, clean reads were produced. At the same time, it was determined what the clean data’s Q20, Q30, GC content, and sequence duplication level were. On clean, high-quality data, all subsequent analyses were performed. During the transcriptome assembly process, all libraries’ and samples’ left files (read1 files) were merged into a single big left.fq file, and all libraries’ and samples’ right files (read2 files) into a single large right.fq file. Using Trinity (Grabherr et al., 2011), transcriptome assembly was completed with min_kmer_cov and all other parameters set to their default settings. Right.fq and left.fq were used as the foundation for this. Gene function was conducted as described before. The alignments of unigenes to the Nt database were carried out using NCBI blast 2.2.28C with an E-value threshold of 1E−5. The program diamond (v0.8.22) was used for the comparison with the Nr (E-value, 1E−5), KOG/COG (E-value, 1E−3), and SwissProt (E-value, 1E−5) databases. Using an E-value limit of 0.01 and hmmscan from HMMER 3.0, Pfam was searched. With an E-value threshold of 1E−6 based on the Nr and Pfam annotations, the GO (Gene Ontology) annotations were carried out in Blast2GO (v2.5) (Wan et al., 2024). Pathway analysis was carried out to ascertain the crucial pathways of DEGs using the databases of the Kyoto Encyclopedia of Gene and Genomes (KEGG) (http://www.genome.jp/kegg). The statistical enrichment of differentially expressed genes in KEGG pathways was evaluated using the KOBAS program. Genes that exhibit substantial levels of differential expression are found in the pathways with an FDR value of 0.05 (Wan et al., 2023).

Five DEGs were chosen for quantitative real-time PCR (qRT-PCR) to validate the transcriptome results. All cDNAs were generated using the Prime Script RT reagent Kit with gDNA Eraser (TaKaRa, Kyoto, Japan), and qRT-PCR experiments were performed using the ABI 7500 Fast Real-Time Detection System and the TaKaRa SYBR Green Mix kit (TaKaRa, Kyoto, Japan). According to Li et al. (2021), the amplification system and process were used. The reference genes were Primer Premier 5.0-designed primers and 18S rRNA. Supplementary Table S1 lists all of the primers utilized. The 2CT approach was used to calculate the relative expression level.

The data were examined using Microsoft Excel 2007. The Statistical Analysis System v. 9.2 was used for significance testing (Duncan’s test) and principal component analysis (PCA). Heatmaps and Venn diagrams were drawn using Tbtools, while other figures were drawn using Origin 2022.

Different concentrations of MeJA affect the accumulation of L-borneol. Overall, different concentrations of MeJA can both promote L-borneol accumulation in the plant’s leaves. Except for the 10 mmol/L MeJA of 120-h treatment, moxa under other concentration treatments, there is accumulation of L-borneol in the leaves of the three leaf positions of Ainaxiang over time, especially when the concentration was 1 mmol/L on 48-h treatment. The treatment effect of 1 mmol/L MeJA was the most effective, resulting in the highest accumulation of L-borneol at 120 h in leaves at different stages of maturity (young, mature, and old leaves), with concentrations of 3.043 mg·g⁻¹ FW, 3.346 mg·g⁻¹ FW, and 2.044 mg·g⁻¹ FW, respectively. These values showed significant differences compared to the control. Figure 1 clearly demonstrates the positive impact of 1 mmol/L MeJA treatment on the accumulation of L-borneol in the leaves of Ainaxiang.

To explore MeJA’s regulation mechanism on terpenoid synthesis in the leaves of Ainaxiang, non-induced and induced conditions were selected. To obtain the transcriptome, Illumina high-throughput next-generation sequencing was performed on the mature leaves of Ainaxiang under different conditions. The clean data screening requirements are tightly regulated in order to assure the caliber of the subsequent analysis. The reads with adapters should be eliminated first, and then, the reads that make up more than 10% but whose base information cannot be identified. The gathered clean reads for the non-reference genome were spliced using Trinity before the sequence was examined. The primary assembly yielded a total of 509,285 transcripts with minimum and maximum lengths of 201 and 23,172, respectively, and most of them (280,106) are between 200 bp and 500 bp, followed by those in the 501–1,000 bp range (82,683 transcripts) ( Supplementary Table S2 ). In addition, 68,679 transcripts were observed in the range of 1,001–2,000 bp, and 77,817 transcripts were in the range of more than 2,000 bp, with average transcript length of 989 bp and N50 of 2,088 bp.

We aligned the 267,178 unigenes to the Nr (NCBI non-redundant protein sequence), Nt (NCBI non-redundant nucleotide sequence), Pfam (protein family), KOG (protein ortholog source clusters), Swiss-Prot (manually annotated and reviewed protein sequence database), and KO (KEGG Ortholog database), achieving a match of 23.47% × 65.36, for example. In Figure 2 , comparing DEG using a volcano blot showed that 593, 224, 612, 2,405, 1,353, and 921 genes were upregulated, and 4, 123, 573, 1,745, 766, and 763 genes were downregulated in the D1_1vsCK, D1_10vsCK, D2_1vsCK, D2_10vsCK, D5_1vsCK, and D5_10vsCK treatments, in that order. We used GO enrichment analysis to anticipate the activities of the DEGs. We found that the DEGs enriched for three primary ontologies: biological process, cellular component, and molecular function. We assigned the 666,068 unigenes to three major GO terms: cellular component (198,296 unigenes across all samples), molecular function (148,057 unigenes), and biological process (319,733 unigenes). Finally, we annotated the 62,720 DEGs of all samples to 25 GO terms, including posttranslational modification, protein turnover, chaperones category, which contained 6,248 unigenes, and translation, which had 2383 unigenes related ( Figure 3 ; Supplementary Table S3 ).

In addition, the KOG database showed that 50,632 unigenes had strong functional prediction and categorization matches ( Figure 4 ; Supplementary Table S4 ). The 25 KOG categories that were looked at were translation, ribosome structure and biogenesis (5963), general function prediction alone (5544), and signal transduction pathways (3005). The largest groups were posttranslational modification, protein turnover, and chaperones (6248). The smallest categories were extracellular structures (46) and cell motility (40). We also examined the unigenes’ KEGG pathways. We found that 25 projected metabolic pathways ( Figure 5 ; Supplementary Table S5 ), grouped into five categories: metabolism, genetic information processing, environmental information processing, cellular activities, and organismal systems, involved 64803 unigenes. The most numerous category in its gene content was metabolism, including carbohydrate metabolism (6244), followed by overview (5105), energy metabolism (4814), amino acid metabolism (3764), and lipid metabolism (3195). As demonstrated in Figure 5 and Supplementary Table S5 , the smallest categories were organismal systems (1772) and environmental information processing (1861).

The expression profiles of genes involved in borneol production were studied to investigate the regulatory mechanisms for the accumulation patterns of diverse borneol in B. balsamifera under different MeJA treatments. B. balsamifera was found to include 4,848 expressed unigenes that encode terpenoid biosynthesis enzymes. At 24 h after 1 mmol/L MeJA treatment, there are 1,077 DEGs, 1,185 DEGs at 48 h, and 2,119 DEGs at 120 h. At each sample point, there were 593, 484, and 612 upregulated genes. They also had 347, 4,150, and 1,684 DEGs at each point under 10 mmol/L MeJA treatment; upregulation genes were 224, 123, and 2,405 at three sample points, respectively. Figure 6 displays the FPKM values and expression data for each unigene. Most genes that code for important enzymes in the borneol backbone pathway (MEP and MVA pathways, KEGG entry ko00900) did have a high transcriptome expression level. However, the deoxy-d-xylulose 5-phosphate synthase (DXS), ispH, geranylgeranyl diphosphate synthase (GGPS), geranyl diphosphate synthase (GPS), and 1-deoxy-d-xylulose-5-phosphate reductoisomerase (DXR) unigenes did not. Under the low concentration of MeJA treatment, DXS (Cluster-28937.97079) exhibits high expression on all treatment days. At the same time, every ispH (Cluster-28937.127564), GGPS (Cluster-28937.10564), and GPS (Cluster-28937.47046) showed upregulation only after 120 h treatment. On the other hand, under 10 mmol/L MeJA treatment, the expression of GGPS (Cluster-28937.10564) and GPS (Cluster-28937.47046) showed upregulation in all treatment days, and DXR (Cluster-28937.158907) showed upregulation at the first 48 h, while DXS (Cluster-28937.97079) exhibited a high expression level after 120 h of treatment. Interestingly, unigenes in the MEP route have greater expression levels than those in the MVA pathway.

The expression analysis related to borneol biosynthesis genes was conducted at both 1 and 10 mmol/L MeJA treatments, as shown in Figure 7 . The gene expression results are consistent with the transcriptomic analysis results. Genes DXR, DXS, and GPS showed high expression over time at 1 mmol/L treatment. At the same time, ispH and GGPS showed high expression with a 10 mmol/L treatment.