A commercial variety of M. arvensis was propagated by planting root segments under natural conditions. Each root segment contained one growth point and was planted in a substrate composed of peat soil and perlite in a 3:1 ratio. Six plants were planted per pot, with pot dimensions of 49 × 21 × 15 cm. The experiment included 3 independent replicates. Approximately 10 weeks after planting, once the plants had developed 5 leaves and entered the vigorous growth phase, the plants were treated with varying concentrations of MJ or SA in order to screen for effective hormone. The MJ concentrations used were 10 μM, 100 μM, 1 mM, 5 mM, and 10 mM, the solution (Aladdin, Shanghai, China) was prepared with 0.2% Tween 20 (Sangon, Shanghai, China). The SA concentrations used were 100 μM and 1 mM, and the solution (Aladdin, Shanghai, China) was prepared with 0.2% Tween 20 (Sangon, Shanghai, China). The control group was treated with a 0.2% Tween 20 aqueous solution. Hormone treatments were applied weekly and sprayed on leaves until there was runoff for a total of 8 treatments. Twenty-four hours after the final treatment, the aboveground parts of the MJ-treated plants were harvested, and measurements were taken for fresh weight, PGT density, EO content, and yield. Similarly, twenty-four hours after the final treatment, the aboveground parts of the SA-treated plants were harvested, and measurements were taken for EO content.

M. arvensis plants were cultured in a temperature-controlled tissue culture room maintained at 25 °C with 70% humidity, under a 16-hour light/8-hour dark cycle with a light density of 150 μmol·m^-2·s^-1. Once the plants developed the first 5 leaves, they were transitioned to continuous light conditions and pre-cultured for one week. The plants were then sprayed with 1 mM MJ solution until runoff. The MJ solution (Aladdin, Shanghai, China) was prepared with 0.2% Tween 20 (Sangon, Shanghai, China). The control group was treated with a 0.2% Tween 20 aqueous solution. The third pair of leaf samples were collected at 0 hours (CK), 4 hours (H4), 8 hours (H8), and 24 hours (H24) after MJ treatment. Three biological replicates were taken for each time point. All samples were immediately frozen in liquid nitrogen and stored at -80°C for further analysis.

The EO of M. arvensis was extracted via hydro-distillation using a Clevenger apparatus. A sample of 200 g of above-ground fresh plants was subjected to hydro-distillation for 50 minutes. EO content was expressed as a percentage of fresh weight (w/w). EO yield was quantified in grams per plant. The composition of the EO was analyzed using GC-MS. An Agilent 7890B gas chromatograph coupled with a 5977A mass spectrometer (Agilent Technologies, Santa Clara, USA) was employed for the GC-MS analysis. The chromatographic separation was performed on a DB-WAX column (30 m × 0.25 mm, 0.25 µm film thickness; Agilent Technologies, Santa Clara, USA). Helium was used as the carrier gas at a flow rate of 1 mL/min. The quadrupole temperature was maintained at 150°C, and the ionization mode was electron ionization (EI+) (Mahmoud and Croteau, 2001).

From our previous preliminary experiments, we observed that the PGT density was the highest in young leaves and the lowest in old leaves. However, despite their higher PGT density, young leaves often contain immature PGTs that have not yet synthesized or stored EO. Consequently, leaves may contain PGTs at various developmental stages and sizes, complicating accurate observation and statistical analysis. The third and fourth pairs of leaves typically have mature and rounded PGTs filled with EO. Therefore, we selected the third pair of leaves of uniform size as representative samples for PGT density analysis. Fluorescence microscopy (Olympus, Tokyo, Japan) was used to capture images of the abaxial side of each leaf, focusing on the top, middle, and bottom sections. A 10x objective lens was used for imaging, and fluorescence was induced using UV light. The quantification of PGT was performed using ImageJ software (version 1.54k, National Institutes of Health, Bethesda, USA) with appropriate image processing and analysis protocols (Yan et al., 2017). This approach ensured precise quantification of trichome density across different leaf sections.

Total RNA from plant leaves was extracted using the Plant RNA Extraction Kit (Takara, Dalian, China) following the manufacturer’s protocol. The purity and concentration of the extracted RNA were assessed with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA). RNA samples with an A260/A280 ratio between 1.8 and 2.0, and an A260/A230 ratio above 2.0 were selected for further analysis. The RNA sequencing was performed by Biomarker Technologies (Qingdao, China). Sequencing libraries were prepared using the NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (New England Biolabs, Ipswich, USA) following the manufacturer’s protocol, with index codes added for sample identification. mRNA was purified from total RNA using poly-T magnetic beads. First-strand and second-strand cDNA synthesis were carried out using M-MuLV Reverse Transcriptase, DNA Polymerase I, and RNase H. After adenylation of 3’ ends, NEBNext Adaptors were ligated. cDNA fragments (~240 bp) were purified using the AMPure XP system (Beckman Coulter, Beverly, USA). The USER enzyme was applied to the adaptor-ligated cDNA before PCR amplification using Phusion High-Fidelity DNA polymerase and specific primers. Finally, PCR products were purified and the library quality was assessed on an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, USA). Clean reads were obtained by processing raw fastq data with in-house scripts to remove adapters, poly-N sequences, and low-quality reads, while calculating Q20, Q30, GC content, and sequence duplication levels for quality control. Clean reads were assembled de novo using Trinity software (version 2.8.4) with default parameters. Functional annotation of the assembled sequences was carried out using several databases, including the NR (NCBI non-redundant protein sequences), Pfam (Protein family), KOG/COG/eggNOG (Clusters of Orthologous Groups of proteins), Swiss-Prot (A manually annotated and reviewed protein sequence database), KEGG (Kyoto Encyclopedia of Genes and Genomes), and GO (Gene Ontology), to provide comprehensive insights into the biological functions of the identified genes (Li et al., 2023).

Differentially expressed genes (DEGs) were identified using DESeq2 software with the threshold false discovery rate (FDR) < 0.05 and |log2 Fold Change| ≥ 1 (Love et al., 2014). The GO classification analysis and KEGG pathway enrichment analysis were performed using the GO database (https://github.com/tanghaibao/Goatools) and KOBAS program (http://kobas.cbi.pku.edu.cn/) with p < 0.05. ClusterProfiler (version 4.12.6, Yulab, Guangzhou, China) was used to visualize the enrichment results using bar graphs (Klopfenstein et al., 2018). Volcano plots were constructed and analyzed using TBtools software (version 2.119, Guangzhou, China), allowing for the identification of significant DEGs by plotting the log2 fold change against the -log10 p-value (Chen et al., 2020). Veen diagram analysis was performed using the VENNY platform (https://bioinfogp.cnb.csic.es/tools/venny/index). Heatmaps were constructed from log2(TPM) values using TBtools (Chen et al., 2020).

TFs were predicted using the Biomarker Cloud platform (http://www.biocloud.net/) (Du et al., 2021). Gene co-expression analysis was conducted via the SRplot platform (http://www.bioinformatics.com.cn) (Tang et al., 2023). DEGs involved in monoterpenoid biosynthesis of M. arvensis, along with key enzyme genes previously reported in monoterpenoid biosynthesis pathways (Lange and Croteau, 1999b; Lange et al., 2011), and differentially expressed transcription factors (DETFs), were used to construct a co-expression trend network diagram ( Supplementary Table S1 ). Eight monoterpenoid biosynthesis enzyme genes ( Supplementary Table S2 ) were mapped to the transcriptome database with a reference genome by local Blast, and the promoter sequence was searched. Promoter binding elements analysis was then performed using the PlantRegMap platform (https://plantregmap.gao-lab.org/) (Du et al., 2021).

DEGs involved in monoterpenoid biosynthesis and DETFs were used to perform weighted gene co-expression network analysis (WGCNA). The WGCNA was conducted using the Biomarker Cloud platform (http://www.biocloud.net/) (Du et al., 2021). Five co-expressed modules were identified by WGCNA. Genes in module 5 were selected to predict and analyze protein–protein interaction (PPI) using the STRING database (https://cn.string-db.org). The interactions were filtered based on a confidence score threshold of 0.7 to ensure high reliability. The resulting PPI network was then imported into Cytoscape (version 3.10.2, Cytoscape Consortium, San Diego, USA). To identify hub genes, the CytoHubba plugin was employed, and the top 15 hub genes were selected based on their MCC scores, indicating their potential key roles in monoterpenoid biosynthesis biological processes under study (Li et al., 2024).

The total RNA of the leaf samples was extracted using the RNAprep Pure Plant Kit (Tiangen, Beijing, China). The quality and concentration of the RNA were assessed by agarose gel electrophoresis (Major Science, Saratoga, USA) and Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA). The cDNA was synthesized using the HiScript II Reverse Transcriptase kit (Vazyme, Nanjing, China). The gene-specific primers were designed with Primer Premier software (version 6, PREMIER Biosoft International, Palo Alto, USA). The specific primers used are shown in Supplementary Table S3 . Quantitative real-time PCR (qRT-PCR) was performed using AceQ qPCR SYBR Green Master Mix kit (Vazyme, Nanjing, China) and a Bio-Rad MiniOpticon Real-Time PCR machine (Bio-Rad, Hercules, USA). The PCR reaction conditions were: pre-denaturation at 95°C for 1 min, followed by denaturation at 95°C for 10 seconds, and annealing at 60°C for 30 seconds, a total of 40 cycles were performed. All the data were normalized using Actin gene as reference, and the gene expression level was calculated using 2^−ΔΔCT (Li et al., 2023).

The values are presented as the mean ± standard deviation (SD) of at least 3 replicates. One-way analysis of variance (ANOVA) was performed using GraphPad Prism (version 8.0, GraphPad Software, San Diego, USA), and statistical significance was set at *P < 0.05, **P < 0.01, ***P < 0.001, and ****P< 0.0001 (Park et al., 2024).

The biosynthesis of metabolites in many medicinal plants is closely linked to defense mechanisms (Bednarek and Osbourn, 2009). MJ and salicylic acid (SA), two key defense-related phytohormones, were applied to mint leaves to evaluate their effects. MJ significantly increased the EO content, while SA showed minimal to no effect ( Supplementary Figure 1 ). Further experiments demonstrated that MJ promoted M. arvensis EO synthesis in a concentration-dependent manner ( Figures 1A–C ). As the MJ concentration increased from 10 μM to 5 mM, a corresponding increase in EO content was observed. At the 5 mM dose of MJ, the EO content peaked, increasing by 71.20% compared to the control group. While high concentrations of MJ slightly reduced plant growth (e.g., a 9.51% decrease in fresh weight under 10 mM treatment), this did not hinder EO production. Notably, under 1 mM MJ treatment, the EO production increased by 58.50% compared to the control group, although the highest EO content was achieved under 5 mM MJ treatment.

PGTs are where mint EO is synthesized and stored (Croteau et al., 2005). The density of mint PGT was assessed using fluorescence microscopy. The results showed that the trends in PGT density were closely aligned with the changes in EO content ( Figures 1D, E ). Specifically, the PGT density peaked under 5 mM MJ treatment, with an increase of 53.69% compared to the control group.

Across all treatments, 28 chemical components were identified in M. arvensis EOs, representing over 98.52% of the total composition ( Table 1 ). In general, the types and relative proportions of the compounds identified remained consistently similar across all treatments. In each group, monoterpenoids were the predominant components. The top 5 compounds were L-menthol (79.07%~85.86%), L-menthone (4.93%~10.89%), neomenthol (1.60%~1.90%), isomenthone (1.20%~1.52%) and pulegone (0.91%~2.66%). These compounds represented over 94.80% of the total volatiles in each treatment. However, subtle variations were observed between different groups, particularly under the 5 mM MJ treatment. L-menthol, the predominant component of M. arvensis EO, serves as a key quality indicator (Adlard, 2010; Wei et al., 2023). In all treatments, L-menthol content exceeded 79.07%. The control group exhibited the highest L-menthol content (85.86%). As MJ concentration increased, L-menthol content decreased slightly. Under 1 mM MJ treatment, L-menthol content was 82.57%, while at a 5 mM MJ dose, it decreased to 79.07%, a significant reduction of 7.90% compared to the control group (p < 0.05). In contrast, under 5 mM MJ treatment, menthone and pulegone content increased significantly by 120.80% (p < 0.05) and 191.72% (p < 0.01), respectively, compared to the control.

Transcriptome sequencing was performed to explore how MJ regulates EO biosynthesis and PGT development in M. arvensis at the transcriptional level. Total RNA was extracted from M. arvensis leaves at 0, 4, 8, and 24 hours after MJ treatment. Twelve cDNA libraries (3 replicates per treatment) were sequenced, yielding 81.56 Gb of clean data. Each library produced no less than 6.41 Gb of clean data. GC content ranged from 46.91% to 47.38% across the 12 transcriptome samples. The percentage of Q30 bases was at least 92.58% ( Table 2 ). Sample replicate correlation analysis ( Supplementary Figure S2A ) showed that the Pearson’s correlation coefficients within the same treatment group were not less than 0.899. Furthermore, the principal component analysis ( Supplementary Figure S2B ) revealed that the values of PC1 and PC2 were 17.14% and 13.23%, respectively. These results indicated that the dataset is robust, with high reproducibility between samples, making it suitable for further analysis.

A total of 55,154 unigenes were annotated using public protein databases (NR, Swiss-Prot, GO, COG, KOG, KEGG, etc.) with a BLAST E-value cutoff of 1.0 × 10^-5 ( Supplementary Table S4 ). To evaluate the effect of MJ on M. arvensis gene expression, DEG analysis was performed. A total of 7,428 DEGs were identified across six comparative groups ( Figure 2 ). In the CK vs H4 comparison, 2,765 DEGs were identified, including 1,317 up-regulated and 1,448 down-regulated genes. In the CK vs H8 comparison, 1,674 DEGs were identified, including 745 up-regulated and 929 down-regulated genes. In the CK vs H24 comparison, 2,443 DEGs were identified, including 1,828 up-regulated and 615 down-regulated genes. In the H4 vs H8 comparison, 2,069 DEGs were identified, including 929 up-regulated and 1,140 down-regulated genes. In the H4 vs H24 comparison, 4,511 DEGs were identified, including 2,971 up-regulated and 1,540 down-regulated genes. In the H8 vs H24 comparison, 3,885 DEGs were identified, including 2,731 up-regulated and 1,154 down-regulated genes ( Figures 2A, B ; Supplementary Figure S3 ). In addition, a total of 4,659 DEGs were identified across the CK vs H4, CK vs H8, and CK vs H24 comparisons. The analysis revealed that 96 genes were consistently up-regulated, while 124 genes consistently down-regulated after MJ treatment ( Figure 2C ). Over time, different genes showed varied expression patterns in response to MJ ( Figure 2D ). Overall, MJ significantly affected gene transcription in M. arvensis.

GO enrichment analysis was performed to explore the potential functions of these DEGs from the CK vs H4, CK vs H8, and CK vs H24 comparisons ( Figure 3 ), categorizing them into 3 major groups: biological process (BP), cellular component (CC), and molecular function (MF). Within BP, the top 4 subcategories were “metabolic process”, “cellular process”, “single-organism process” and “biological regulation”. In CC, the top 3 subcategories were “membrane”, “cell” and “cell part”. In MF, most DEGs clustered in “binding”, followed by “catalytic activity” ( Figure 3A ). Notably, 1,569 DEGs were annotated in the “metabolic process” subcategory, making it the most abundant category in the GO enrichment analysis. Additionally, 409, 203, 45, 41, and 20 genes were annotated in the subcategories of “response to stimulus”, “detoxification”, “immune system process”, “signal transducer activity”, and “antioxidant activity”, respectively ( Figure 3A ; Supplementary Table S5 ). These subcategories have been shown to play key roles in plant resistance to stress.

KEGG pathway analysis showed that 130 pathways were significantly enriched under MJ treatment. DEGs were predominantly enriched in the plant-pathogen interaction (240 genes), MAPK signaling (179 genes), phenylpropanoid biosynthesis (134 genes), and plant hormone signal transduction pathways (126 genes). Notably, 28 genes were enriched in the monoterpenoid biosynthesis pathway ( Figure 3B ; Supplementary Table S6 ). Additionally, 12 genes (1.17%), 9 genes (1.53%), and 19 genes (2.13%) were enriched in this pathway in the CK vs H4, CK vs H8, and CK vs H24 comparisons, respectively ( Supplementary Figure S4 ), indicating that H24 induced a higher number of genes involved in monoterpenoid biosynthesis. For example, in the CK vs H24 comparison, 3 genes encoding cytochrome P450 enzyme (CYP450) (TRINITY_DN14738_c1_g1, TRINITY_DN32509_c0_g1, and TRINITY_DN46431_c0_g1), 2 genes encoding IPR (TRINITY_DN60151_c0_g1 and TRINITY_DN90669_c0_g1), 1 gene encoding LS (TRINITY_DN5134_c0_g1), and 1 gene encoding MD (TRINITY_DN28821_c1_g1 encoded) were identified.

In Arabidopsis thaliana, the F-box protein coronatine insensitive 1 (COI1) acts as the receptor for JA. Perception of JA-Ile by the SCF^COI1 complex triggers the degradation of JAZ proteins through the 26S proteasome. This process activates downstream TFs involved in JA responses, such as Octadecanoid responsive 3 (ORA3), a member of the APETALA2/ethylene responsive factor (AP2/ERF) family (Montiel et al., 2011; Yu et al., 2021b). In this study, 10 DEGs were identified as key factors of JA signaling in M. arvensis ( Figure 4A ). Among them, 1 JAR6 gene (TRINITY_DN3071_c0_g1) and 4 JAZ9 genes (TRINITY_DN16183_c0_g1, TRINITY_DN19929_c0_g1, TRINITY_DN11251_c0_g1 and TRINITY_DN834_c0_g2) were significantly up-regulated at 4 hours after MJ treatment. Two MYC2 genes (TRINITY_DN4107_c0_g2 and TRINITY_DN16914_c0_g1) were significantly up-regulated at 24 hours after MJ treatment ( Figure 4A ; Supplementary Table S1 ).

In M. arvensis, monoterpenoids are the primary components of its EOs. KEGG enrichment results indicated that “Monoterpenoid biosynthesis” was one of the most significantly enriched pathways under MJ treatment. Twenty-eight DEGs associated with monoterpenoid biosynthesis were identified in response to MJ treatment. Notably, GPPSs were significantly up-regulated at 4 and 8 hours after MJ treatment, 4 LS genes were up-regulated both at 4, 8 and 24 hours, and 12 IPR genes were significantly up-regulated at 24 hours after MJ treatment. CYP450, catalyze the decoration of terpenoid basic skeletons and thereby contribute significantly to their structural diversity (Weitzel and Simonsen, 2015). In this study, 4 CYP450 genes were identified as significantly expressed ( Figure 4B ; Supplementary Table S1 ). L-menthol is the predominant component of M. arvensis EO, which is synthesized through a series of enzymatic reactions. The enzymatic catalysis mechanism of L-menthol synthesis has been extensively studied in 3 Mentha varieties: Mentha piperita, Mentha spicata, and Mentha haplocalyx (Ahkami et al., 2015). Using reference genes from these varieties as queries, 11 orthologous genes were identified in M. arvensis ( Supplementary Table S1 ), including the previously reported LS and MD genes (Akhtar et al., 2017; Wang et al., 2013).

TFs regulate plant growth, development, and secondary metabolite synthesis (Yang et al., 2012). In this study, 260 DETFs belonging to 22 families were identified after MJ treatment for 4, 8, and 24 hours ( Figure 5A ; Supplementary Table S1 ). The top 8 families with the highest number of DEGs were AP2/ERF (42 genes), MYB (including MYB-related, 38 genes), WRKY (33 genes), NAC (22 genes), HSF (heat shock transcription factors, 16 genes), bHLH (14 genes), C2H2 (C2H2 zinc-finger protein, 14 genes), and GRAS (GAI-RGA-and-SCR, 14 genes). These TFs responded positively to MJ treatment. To further explore the regulatory relationship between TFs and monoterpenoid biosynthesis, a co-expression pattern analysis of DETFs and monoterpenoid biosynthetic genes was performed ( Figure 5B ). The analysis identified 7 gene clusters. Each cluster exhibited a unique expression pattern, indicating a close association between TFs and monoterpenoid biosynthetic enzyme genes. Most of the highly expressed monoterpenoid biosynthesis genes (Log2FC > 3) were grouped into cluster 6. Therefore, further analysis focused on this cluster. In cluster 6, the top 3 TF families with the highest number of DEGs were AP2/ERF (17 genes), WRKY (16 genes), and MYB (10 genes) ( Figure 5B ; Supplementary Table S2 ). In this study, 8 monoterpenoid biosynthesis genes, including DXS, DXR, GPPS-L, LS, IPR, PR, and MD, which are either highly expressed genes (Log2FC > 3) or key genes in this process as previously reported (Lange and Croteau, 1999a; Lange et al., 2011; Mahmoud et al., 2004) were selected for promoter analysis. These genes contained the highest number of AP2/ERF binding sites (2,763), followed by bHLH (266), WRKY (251), and MYB (197) ( Figure 5C ; Supplementary Table S2 ).

WGCNA was employed to construct a gene co-expression network, aiming to identify DETFs involved in MJ-induced monoterpenoid biosynthesis. Hierarchical clustering identified 5 co-expressed modules ( Figure 6A ). Most of the high expressed monoterpenoid biosynthesis genes (10 out of 11, Log2FC > 3) peaked at 24 hours after MJ treatment and were classified in module 5 ( Figure 6B ; Supplementary Table S7 ). KEGG analysis was performed for this module. Among these, 18 unigenes were annotated in the plant hormone signal transduction pathway (35.29%), 16 unigenes in the plant pathogen interaction pathway (31.37%), and 11 unigenes in the monoterpenoid biosynthesis pathway (21.57%) ( Supplementary Figure S4 ). Additionally, module 5 was highly associated with the H24 group (R = 0.99, p = 0.006). These results closely mirrored the analysis of the monoterpenoid biosynthesis pathway in Figure 4 , prompting a focus on module 5 ( Figure 6B ; Supplementary Table S7 ). Furthermore, all genes of module 5 were used to conduct PPI analysis, the results were visualized using Cytoscape, and hub genes were identified using CytoHubba ( Figure 6C ). A venn diagram analysis was used to identify common genes between the top 15 hub genes from the CytoHabba analysis and cluster 6 genes from the co-expression analysis (Log2FC > 1, FPKM > 10) ( Figure 5B ), which were considered as candidate genes. The results revealed 5 common genes, including 4 AP2/EFR genes—TRINITY_DN8517_c0_g1 (ERF108), TRINITY_DN4459_c0_g1 (ERF1B), TRINITY_DN11107_c0_g1 (DREB1D), TRINITY_DN688_c0_g1 (DREB1C), and 1 WRKY gene TRINITY_DN1548_c0_g1 (WRKY33), which were significantly up-regulated (Log2FC > 2) at 24 hours after MJ treatment ( Figures 6D, E ).

A qRT-PCR assay with independent samples from the control and MJ treatment group (H24) was conducted to verify the expression changes of several key genes involved in monoterpenoid biosynthesis. In total, 10 genes, including 6 from the monoterpenoid biosynthesis pathway, and 4 TFs were selected to confirm the RNA-seq data. The expression levels of these selected genes, as determined by qRT-PCR, were generally consistent with the RNA-seq results ( Figure 7 ).