Transcriptome analysis of transcription factors and enzymes involved in monoterpenoid biosynthesis in different chemotypes of Mentha haplocalyx Briq

Plant materials

M. haplocalyx samples (l-menthol chemotype, pulegone chemotype, and carvone chemotype) from the Medicinal Botanical Garden of Jiangsu University, in Zhenjiang City, Jiangsu Province, with geographic coordinates of N32°12′5″, E119°30′3″, were collected approximately every 10 days from May 10, 2020, from 12:00 to 14:00 for 5 months. Each sample had three biological replicates. The above-ground parts were dried at room temperature for the determination of essential oil content. The top leaves were collected, immediately wrapped in aluminum foil and frozen in liquid nitrogen, and then stored at −80 °C for transcriptome sequencing and chemotype variation research. Table S1 shows the samples and local temperature.

Determination of characteristic terpenoids in essential oils

The preparation of essential oil followed the method outlined in Chinese Pharmacopoeia (2020 edition) (Chinese Pharmacopoeia Commission, 2020). The contents of the essential oils were determined using the method used by our group in a previous study using gas chromatography (GC) (Yang, 2020; Yang et al., 2021). GC conditions were as follows: the DB-1 column was 60 m × 0.25 mm × 0.25 μm in size and the inlet and detector temperature were both 250 °C. The split ratio was 40:1, the carrier gas was N₂ (99.999%), and the carrier gas flow rate was 1 mL/min. The heating program for l-menthol chemotype and pulegone chemotype was as follows: the initial temperature was 50 °C for 28 min, increased to 90 °C at 10 °C/min, increased again to 150 °C at 3 °C/min, increased to the final 230 °C at 10 °C/min, and held for 5 min. The heating program for carvone chemotype was as follows: the initial temperature was 60 °C for 4 min, increased to 80 °C at 2 °C/min, increased to 150 °C at 3 °C/min, then increased to 230 °C at 5 °C/min, and held for 5 min. All determinations were performed in triplicate.

cDNA library construction and transcriptome sequencing

The stable l-menthol chemotype, pulegone chemotype, and carvone chemotype samples collected on August 5^th, 2020 and pulegone chemotype sample collected on May 20^th, 2020 were selected for transcriptome sequencing. There were three biological replicates for each type. Independent total RNA extracts were prepared from the 12 samples and equal amounts of RNA from each sample were used to generate RNA-Seq libraries. The total RNA of M. haplocalyx chemotypes was extracted using Trizol reagent (Sangon Biotech, Shanghai, China). The quality was assessed using 1% agarose gel electrophoresis and a Qubit 2.0 RNA kit (Life Technologies, Carlsad, CA, USA). Table S2 shows the RIN value. A total amount of 2 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using Hieff NGS™ MaxUp Dual-mode mRNA Library Prep Kit for Illumina® (Yeasen Biotechnology, Shanghai, China) following manufacturer’s recommendations. First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-) (Yeasen Biotechnology, Shanghai, China). Second-strand cDNA synthesis was subsequently performed using DNA polymerase and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, adaptor was ligated to prepare for library generation. In order to select cDNA fragments of a preferred 150–200 bp length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Then, 3 μL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. PCR was performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. Finally, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The libraries were then quantified and pooled. Paired-end sequencing of the library was performed on the Illumina Novaseq 6000 (Illumina, San Diego, CA, USA).

Transcriptome assembly and annotation

The original reads obtained by RNA sequencing were filtered with Trimomatic (version 0.36) to remove low-quality reads and the remaining clean reads were reassembled into transcripts by Trinity (version 2.0.6) with default settings to obtain unigenes for further analysis. Transcripts with a minimum length of 200 bp were clustered to minimize redundancy. All data generated were saved in the National Center for Biotechnology Information (NCBI) and can be accessed in the Short Read Archive (SRA) sequence database (accession number PRJNA795820). An online implementation of RNASeqPower is available at https://rodrigo-arcoverde.shinyapps.io/rnaseq_power_calc/. Unigenes were blasted against NCBI Nr (NCBI non-redundant protein database), SwissProt (A manually annotated and reviewed protein sequence database), TrEMBL (A supplement to the SwissProt database), CDD (Conserved Domain Database), PFAM (Protein family), and KOG (Eukaryotic Orthologous Groups) databases (E-value < 1e−5).

Analysis of differentially expressed genes

In RNA-seq analysis, we introduced Transcripts Per Million (TPM) to measure the proportion of a certain transcript in the RNA pool. We used Salmon (version 0.8.2) to calculate unigene readings and TPM. DESeq (1.26.0) was used to analyze differential expression. In order to obtain the significant differential gene, the screening conditions were set as follows: Q value < 0.05 and difference multiple |FoldChange|>2. We have chosen the groups M8 vs P8, M8 vs C8, and M8 vs P8 for analysis.

Quantitative real-time PCR

The qRT-PCR method was used to verify the expression of selected genes in M. haplocalyx. The BeyoRT^TM II kit was used for first-strand cDNA synthesis (Beyotime Biotechnology, Shanghai, China) from total RNA. AceQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Nanjing, China) was used to perform the qRT-PCR reaction. qRT-PCR was conducted using the Light Cycler 96 Real-Time PCR system (Roche, Basel, Switzerland) with a 20 µL reaction mixture. The mixture consisted of 10 µL 2 × AceQ Universal SYBR qPCR Master Mix, 2 µL primer mix, 2 µL diluted cDNA, and 6 µL ddH2O. The qRT-PCR program comprised a preincubation step at 95 °C for 600 s followed by PCR: 40 cycles at 95 °C for 10 s and 60 °C for 30 s followed by a melting cycle: 95 °C for 15 s, 60 °C for 60 s and 95 °C for 15 s. Primer 3.0 online software (http://primer3.ut.ee/) was used to design primers. Expression levels of differentially expressed genes (DEGs) were determined by the method of 2^−ΔΔCt (Liu et al., 2020). We selected the β-actin gene as the reference gene (Qi et al., 2018). All results are representative of three independent experiments. Table S3 lists the qRT-PCR primers.

Data analysis

The data are expressed as mean ± standard deviation (SD) and differences were analyzed by one-way ANOVA. The visualization of correlation analysis data was performed using Cytoscape (version 3.9.1) software.

Dynamic changes of characteristic terpenoids in essential oil of M. haplocalyx chemotypes

As shown in Fig. 2, the content of characteristic monoterpene components of different chemotypes was low and unstable until July. With the growth of M. haplocalyx, the chemotypes showed a stable trend and l-menthol, pulegone, and carvone content levels were all high in different chemotypes in August. However, in May, the menthone content was higher in pulegone chemotype. During the growing period, (−)-menthol was the main component in l-menthol chemotype, accounting for 44.20–51.71% of the total essential oil (Fig. 2A). During the harvest period (July to October), pulegone was the main component in pulegone chemotype, accounting for 32.41–54.10% of the total essential oil (Fig. 2B). Carvone was the main component in carvone chemotype during the growing period, accounting for 65.53–79.93% of the total essential oil (Fig. 2C).

Transcriptome sequencing and analysis

We selected stable l-menthol chemotype (M8), pulegone chemotype (P8), and carvone chemotype (C8) samples collected in August 2020 and pulegone chemotype (P5) samples collected in May 2020 to construct 12 cDNA libraries. Transcriptome sequencing depth was 6G. The result of RNA-Seq Power analysis was 0.84. In the transcriptome data, the Q20 and Q30 values of each sample were greater than 98.41% and 94.53%, respectively. The GC content of the 12 samples was 49.14–55.01%. In addition, the mapped reading of each sample was not less than 92.48%, indicating high-quality RNA-seq data (Table 1).

In total, 254,942 unigenes were generated. The unigenes obtained showed a total length of 128,420,176, a mean length of 503.72 bp, and a mean N50 of 621 bp (Table S4). There were 103,991 unigenes annotated in the Nr database, accounting for 40.79% of total genes; 88,978 (34.90%) unigenes annotated in the SwissProt database; 52,541 (20.61%) unigenes annotated in the KOG database; and 9,252 (3.63%) unigenes annotated in the KEGG database. In total, there were 121,796 (47.77%) unigenes that could be annotated in at least one database (Table 2).

Differentially expressed genes and enrichment analysis

From the statistical results, in the comparisons of M8 vs P8, M8 vs P5, and M8 vs C8, 9,676 (5,747 upregulated, 3,929 downregulated), 11,551 (7,663 upregulated, 3,888 downregulated), and 3,432 (2,026 upregulated, 1,406 downregulated) DEGs were obtained, respectively (Fig. S1). The numbers of upregulated and downregulated genes in the three groups (M8 vs P8, M8 vs P5, M8 vs C8) were 1,139 and 670, respectively (Fig. 3).

Through GO enrichment analysis to identify the biological functions of up- and down-regulated DEGs in different comparison combinations, we found that the significant enrichment of upregulated and downregulated DEGs was divided into three categories: biological processes, cell composition, and molecular function (Fig. 4; Table S5). The top 30 pathways with the highest degree of enrichment from the KEGG enrichment analysis are shown in Fig. 5 (Table S6). Terpenoid backbone biosynthesis was significantly enriched in the three comparison groups; thus, this pathway may have a crucial function in the variation of different chemotypes. For the DEG intergroup comparison scheme used in Figs. 4 and 5, we chose the groups M8 vs P8, M8 vs C8, M8 vs P8.

Identification and analysis of DEGs in monoterpenoid biosynthetic pathways of different chemotypes

Ten unigenes were selected for qRT-PCR to verify the reliability of RNA-seq data. Fig. S2 shows that the qRT-PCR expression patterns of most selected unigenes closely matched the TPM results of RNA-seq data, which indicated that the data could be used for subsequent analyses. DEGs of key enzymes were analyzed using the transcriptome annotation results and the main monoterpenoid biosynthetic of the KEGG metabolic pathway. Significant differences were found for genes GPPS, LS, L3OH, iSPD, PR, MD, NMD, and L6OH. Statistics are shown in Table S7.

All relative expression analyses of the selected unigenes were performed based on M8, the l-menthol chemotype. In the comparison groups of different chemotypes, we observed that the relative expression of enzyme genes in the monoterpenoid biosynthesis pathway was related to the main components of different chemotypes of essential oil (Fig. 6). In pulegone chemotypes, the expression level of PR (TRINITY_DN81154_c1_g1 and TRINITY_DN81915_c0_g3) was significantly lower than for other chemotypes (P < 0.05). However, the relative expression level of MD (TRINITY_DN67710_c6_g2 and TRINITY_DN89117_c1_g2) in l-menthol chemotypes was significantly greater than that of other chemotypes (P < 0.05). Highly expressed MD may enable the biosynthesis of monoterpenes to eventually generate l-menthol. L6OH (TRINITY_DN74923_c4_g1) expression in carvone chemotypes was greater compared to that in other chemotypes. The major component of P8 was pulegone, the major component of P5 was menthone, and the major component of M8 was (−)-menthol. Both menthone and menthol are downstream of PR in the monoterpenoid biosynthesis pathway. Fig. 6 shows that the relative expression of PR in P8 is lower than that in M8, inferring that PR maintains an important role in M. haplocalyx. Similarly, MD, L6OH, and L3OH are also enzymes with important positions in this pathway in M. haplocalyx. It could be argued that the abundance of enzyme expression might affect the variation of chemotypes.

Identification and analysis of the differentially expressed TFs of different chemotypes

In total, 2,599 TFs from 66 families were identified in M. haplocalyx RNA-seq data. The C2H2, MYB, C3H, AP2/ERF, bHLH, WRKY, GRAS, MYB, bZIP, and NAC TF families were the 10 largest families (Fig. S3). The differential TFs identified in the transcriptome data included 113 TFs from 34 families, among which GRAS, bHLH, C2H2, ERF, Dof, MYB, TCP, B3, WRKY, C3H, and Trihelix contained more than five DEGs.

Interaction network analysis of genes related to monoterpenoid biosynthesis of different chemotypes

To further investigate the regulatory mechanism of monoterpenoid biosynthetic pathways of different chemotypes, we used Cytoscape to visualize the correlation between the key enzyme genes of monoterpenoid biosynthetic pathways and the expression profiles of TF genes (Fig. 7; Table S8). According to a past study (Dong et al., 2020), the main TFs involved in the regulation of terpenoid synthesis include six TF families: AP2/ERF, bHLH, MYB, NAC, bZIP, and WRKY. It is worth noting that in the comparison of M8 vs. P8, bHLH, bZIP, ERF, and MYB showed very high correlations with MD and PR (r > 0.95) (Fig. 7A). The expression of eight unigenes of these TFs families was significantly different (P < 0.05, Fig. 7D). In the comparison of M8 vs. P5, bHLH, bZIP, AP2, MYB, and WRKY showed very high correlations with MD and PR (r > 0.95) (Fig. 7B). The expression of seven unigenes of the bHLH, bZIP, MYB, and WRKY families was significantly different (P < 0.05, Fig. 7D). We speculated that these TFs might affect the essential oil components in different chemotypes by regulating the expression of MD and PR. The expression patterns of these key enzymes in turn regulate the variation of l-menthol and pulegone chemotypes. In the comparison of M8 vs. C8, ERF showed a high correlation with L3OH and L6OH (r > 0.85) (Fig. 7C). ERF (TRINITY_DN78234_c2_g1) was significantly differentially expressed (P < 0.05, Fig. 7D). Therefore, these TFs may regulate the transcription of L3OH and L6OH and thus their chemotypes. The relative expression trends of related TFs in different chemotypes were verified, as shown in Fig. 7D. All relative expression analyses of the selected unigenes were performed based on M8, the l-menthol chemotype. Although ERF (TRINITY_DN78234_c2_g1), bHLH (TRINITY_DN70782_c2_g1 and TRINITY_DN85127_c2_g1), bZIP (TRINITY_DN84983_c2_g1), MYB-related (TRINITY_DN72105_c0_g1) were all significantly different compared to M8 (P < 0.05), some of them were not highly correlated with key enzymes (Table S8).