Systematic Identification of Suitable Reference Genes for Quantitative Real-Time PCR Analysis in Melissa officinalis L

Melissa officinalis L. is well known for its lemon-scented aroma and various pharmacological properties. Despite these valuable properties, the genes involved in the biosynthetic pathways in M. officinalis are not yet well-explored when compared to other members of the mint family. For that, gene expression studies using quantitative real-time PCR (qRT-PCR) are an excellent tool. Although qRT-PCR can provide accurate results, its accuracy is highly reliant on the expression and stability of the reference gene used for normalization. Hence, selecting a suitable experiment-specific reference gene is very crucial to obtain accurate results. However, to date, there are no reports for experiment-specific reference genes in M. officinalis. Therefore, in the current study, ten commonly used reference genes were assessed for their suitability as optimal reference genes in M. officinalis under various abiotic stress conditions and different plant organs. The candidate genes were ranked based on BestKeeper, comparative ΔCt, geNorm, NormFinder, and RefFinder. Based on the results, we recommend the combination of EF-1α and GAPDH as the best reference genes to normalize gene expression studies in M. officinalis. On the contrary, HLH71 was identified as the least-performing gene. Thereafter, the reliability of the optimal gene combination was assessed by evaluating the relative gene expression of the phenylalanine ammonia lyase (PAL) gene under two elicitor treatments (gibberellic acid and jasmonic acid). PAL is a crucial gene involved directly or indirectly in the production of various economically important secondary metabolites in plants. Suitable reference genes for each experimental condition are also discussed. The findings of the current study form a basis for current and future gene expression studies in M. officinalis and other related species.


Introduction
Melissa officinalis L., commonly known as the lemon balm is a medicinal herb belonging to the mint family (Lamiaceae). It is well known for its characteristic lemon-scented aroma. It has applications in several industries such as agriculture, sanitary, food, and beverages [1,2]. Traditionally, it has been used to remedy gastrointestinal disorders and several other conditions such as sore throat, cough, insomnia, vertigo, nausea, rheumatoid arthritis, and neurological disorders [1,2]. Additionally, several recent studies have also demonstrated M. officinalis to possess hypolipidemic, antioxidant, antimicrobial, anxiolytic, anticancer, and anti-inflammatory effects [2][3][4]. These pharmacological properties of medicinal and aromatic plants are mainly attributed to the secondary metabolites and essential conditions are shown in Figure 1. In the case of all the designed primers, a single band (in 1.8% agarose gel) and a single peak (qRT-PCR amplification) were obtained (Figures 2  and 3, respectively). The Ct values of the 10 nominees exhibited a wide-ranging value in all samples. Irrespective of the experimental conditions, F1-ATPase showed the lowest Ct value, indicating high expression of the gene.

Expression Stability of Candidate Reference Genes by Comparative ∆Ct and Bestkeeper
The comparative ∆Ct and BestKeeper were used to calculate the stability of the candidate reference genes based on the average standard deviation and the standard deviation (+/-crossing point) values, respectively. The expression stability of candidate genes analyzed by ∆Ct and BestKeeper is shown in Table 2. According to the comparative ∆Ct analysis, GAPDH can be considered the most stable reference gene for plant organs and salt stress. However, in these stresses, BestKeeper recommends using GST. In the cases of heat and osmotic stress, comparative ∆Ct analysis is suggested using 28S rRNA whereas according to the BestKeeper analysis, HLH71 is the best candidate. For cold stress, comparative ∆Ct and BestKeeper analyses recommend using EF-1α and 60S rRNA, respectively. In the cases of all stress combined, 28S rRNA and GST were indorsed by comparative ∆Ct and BestKeeper analyses, respectively. For experiments with in vitro and all conditions combined, comparative ∆Ct analysis proposes using β-tubulin and 28S rRNA, correspondingly. According to BestKeeper analysis, for experiments under in vitro conditions and all conditions combined, HLH71 and F1-ATPase (respectively) were highly recommended. Expression levels of the candidate genes is represented by three main colors (green, black, and red). The red color indicates a high Ct value or low expression and the green color indicate a low Ct value, which corresponds to high expression. The black color represents neither high nor low Ct values (moderate values).

Expression Stability of Candidate Reference Genes by NormFinder
NormFinder software ranked the candidate genes based on their stability values. This software first assessed the intra-and inter-group variations, which are then combined into a direct variation value, which is also known as stability value. The gene with the least variation is finally ranked as the best by the software. The expression stability of candidate genes analyzed by NormFinder is shown in Table 3. According to the NormFinder analysis, GAPDH is the best reference gene for all the experimental conditions tested in this manuscript, except for experiments under osmotic stress and in vitro conditions. For osmotic stress, the software suggests using 28S rRNA whereas under in vitro conditions, using EF-1α is highly recommended.

Expression Stability of Candidate Reference Genes by GeNorm
In geNorm analysis, the software evaluates the gene expression values according to their respective M values. The findings of geNorm analyses are shown in Figure 4. According to our analysis, the most stable reference gene varied across different experiments. In the case of experiments with plant organs, geNorm recommends using a combination of GAPDH and HLH71, whereas, for experiments under heat stress, a combination of ATP-synthase and β-tubulin is recommended. In the case of cold, salt, and osmotic stress, geNorm analyses suggest using a combination of EF-1α and 28S rRNA, GAPDH and βtubulin and ATP-synthase and 28S rRNA, respectively. When all stresses are considered, this software recommends using a combination of GAPDH and 28S rRNA. For gene expression experiments with in vitro M. officinalis plants, geNorm analyses recommend using a combination of ATP-synthase and HLH71. When all experimental conditions are taken together, geNorm analyses using a combination of GAPDH and ATP-synthase. Additionally, we also identified the optimal number of genes required for each experimental condition. The cut-off M-value for the optimal number of required genes is 1.5. In the present study, the pairwise variation value of all the datasets (V n /V n+1 ) was found to be higher than the recommended value, i.e., 0.15 ( Figure 5). Nevertheless, these values are only indicative, and using a higher number of reference genes would increase experimental complexity and instability, eventually leading to ambiguous results. Considering that the V n /V n+1 value is indicative and not a mandatory criterion, one or a combination of two reference genes could be adequate for the normalization of the qRT-PCR data in M. officinalis.

Comprehensive Ranking of the Candidate Reference Genes by RefFinder
RefFinder web-based tool combines geNorm, NormFinder, BestKeeper, and Ct values, and provides a comprehensive ranking. According to the RefFinder, GAPDH can be deemed as the most stable reference gene for experiments with plant organs and under salt stress. The expression stability of candidate genes analyzed by RefFinder is shown in Table 4. In the case of experiments under heat stress, osmotic stress, and all stresses combined, 28S rRNA is the most stable reference gene. For the gene expression experiments under cold stress and for in vitro conditions, RefFinder endorses using GST and β-tubulin, respectively. While considering all the experiments together, RefFinder recommends using 28S rRNA.

Identification of the Most Suitable Reference Gene and Validation Experiments
Based on comprehensive analysis by RefFinder, a combination of EF-1α and GAPDH was identified as the most appropriate reference gene for gene expression studies in Melissa officinalis. On the contrary, HLH71 as a reference gene for gene expression studies in lemon balm is not recommended. Suitable genes for each experimental condition were also selected. For experiments under plant organs and heat stress, ranking suggests using a combination of GAPDH and β-tubulin. In the case of experiments under cold stress, our analysis recommends using a combination of EF-1α and 28S rRNA, whereas, for experiments under salt stress, we endorse using a combination of GAPDH and βtubulin. For experiments under osmotic stress and all stresses combined, we suggest using a combination of ATP-synthase and 28S rRNA and GAPDH and 28S rRNA, respectively. In the case of experiments with in vitro organs, we recommend using EF-1α and β-tubulin. Further, to substantiate the reliability of the candidate gene, the relative expression of the PAL gene under two different elicitor stresses was evaluated. Both the best and the least stable candidate genes were used for normalization. Irrespective of the reference gene type and elicitor type (gibberellic acid and jasmonic acid), we found elevated relative gene expression levels of PAL ( Figure 6). In the case of both elicitors, when normalized with the best reference gene, the PAL gene shows~2X overexpression compared to the untreated control. However, when normalized with HLH71, the PAL gene showed~10X and~15X overexpression under gibberellic acid and jasmonic acid stress, respectively.  stress (B,D). Relative gene expression before elicitor treatment (Control) and after treatment (the most stable and the least stable gene) were compared, and normalization was performed with EF-1α|GAPDH (the most stable gene) and HLH71 (the least stable gene). The relative gene expression between the control and elicitor-stressed samples was compared by Student's t-test, "*" on the vertical bars indicating a significant difference at p < 0.05.

Discussion
Gene expression studies using qRT-PCR have been used across numerous research studies to unravel the molecular mechanisms behind the biosynthetic pathways of various medicinal and aromatic plants [6,[28][29][30][31]. With substantial progress, many genes and gene families have been identified that are directly or indirectly responsible for the synthesis of economically important secondary metabolites [31][32][33]. To identify these genes, selecting suitable reference genes for qRT-PCR analyses is crucial to avoid misleading results. The expression and stability of any given reference gene could vary greatly under different experimental conditions [17][18][19]. Consequently, several studies have identified experimentspecific reference genes in various plant species [18,19]. Despite this progress across numerous plant species, it is surprising that a suitable reference gene in M. officinalis has never been established, in any experimental condition. To date, gene expression studies in M. officinalis have mostly (virtually all the studies) utilized different random reference genes, presuming the selected gene has a stable expression. Hence, the current study is the only study to systematically identify the most suitable reference genes in lemon balm.
In this study, we assessed the stability of ten candidate genes under various experimental conditions to find the best reference gene in M. officinalis. For this purpose, five commonly used software and algorithms for reference gene selection were used. Based on the pairwise variation analysis, we found that a combination of two reference genes will be sufficient for normalization. The geNorm and the NormFinder are specialized software, which can predict the best combination based on expression stability values. According to geNorm, under all the combined stress conditions, the best pair is GAPDH and ATP-synthase, whereas, under similar conditions, NormFinder recommends using a combination of EF-1α and GAPDH. However, the comprehensive ranking software, RefFinder recommends not using ATP-synthase in lemon balm. Finally, based on our comprehensive analyses, a combination of EF-1α and GAPDH was found to be the best reference gene pair whereas HLH71 was found to be the least-performing gene. EF-1α plays important roles during several important processes required for cell growth and proliferation, apart from its role in translation, elongation, and the nuclear export of proteins [34]. Similarly, GAPDH also plays important roles in several non-metabolic processes, such as transcription activation, apoptosis, and axoplasmic transport [35]. Hence, considering their important roles, these genes are expected to be ubiquitously expressed. Moreover, our results are consistent with a previous study where EF-1α and GAPDH along with actin were found to be the most suitable reference gene combination in two thyme species [36]. Similarly, GAPDH was also found to be a suitable reference gene in Rosmarinus officinalis L. and Salvia hispanica L. [31,37]. In Isodon rubescens (another member of the lamiaceae family), using a combination of GAPDH, 18S, and eIF was recommended [38]. However, these results are not consistent across the lamiaceae and other families, for instance, EF1 and EF-1α have been reported to be the least-performing genes in some members of the lamiaceae family [39,40]. In another study, a combination of F1-ATPase, ATP synthase, and ACCase was ranked as the best reference gene for gene expression studies in in vitro-produced rosemary [41]. For Populus euphratica, the optimal reference genes varied for different treatments (e.g., RPL17 was selected for ABA and EF-1α was selected for cold) [42]. In blueberry, GAPDH, ATP1, NADH, and COX2 genes showed high stability [43]. These studies evidently suggest that no reference gene can be considered a universal gene.
Further, the best and least gene was used to normalize the gene expression of the PAL gene under two elicitor treatments (gibberellic acid and jasmonic acid) to affirm its reliability. It is well established that plant hormone treatment such as gibberellic acid and jasmonic acid, stimulates secondary metabolites production by influencing the genes involved in phenylpropanoid pathway of the plant [6]. Correspondingly, in the current study, the expression of the PAL gene was elevated by the hormone treatments. However, in the case of normalization with the least stable gene, the expression levels of PAL were abnormally high, which clearly indicates the impact of using an unsuitable reference gene. Previously, a similar type of result was obtained in a study under herbicide stress by Sen et al. [12], where an herbicide target gene showed abnormal expression when treated with herbicide. Similarly, in our case, treatment with the elicitors is expected to increase the expression of the PAL gene, but in the case of normalization with the least stable reference genes, the expression pattern is misleading.
Overall, the reference genes recommended in this study can be used for gene expression studies under various developmental stages and diverse abiotic stresses in Melissa officinalis and other related species. We recommend using our reference genes to avoid any misleading results in qRT-PCR experiments with lemon balm and other related species.

Plant Material Acquisition, Stress, and Elicitor Treatment
Seed samples of M. officinalis were collected from the well-maintained botanical garden of the Faculty of Tropical AgriSciences (FTZ), Czech University of Life Sciences, Prague (50 • 07 52.9" N, 14 • 22 14.7" E). No prior specific permissions were necessary to collect the seed samples. The seeds were then sown in plastic pots (9 × 9 cm 2 ) containing garden soil (Agro profi RS1, 's-Gravenzande, Netherlands) and perlite in a 3:1 ratio until germinated. Thereafter, the seedlings were maintained in greenhouse conditions with an average temperature of 25 • C and relative humidity between 60 and 70% for 30 days before subjecting to treatments. Five pots per treatment and one plant in each pot were maintained for 24 h. The treatments included: heat stress (35 • C), cold stress (4 • C), salt stress (NaCl, 200 mM), drought stress (PEG 6000, 200 mM), jasmonic acid (1 mg/L), gibberellic acid (1 mg/L) and control (plants without any treatment). The heat and cold stresses were carried out in a growth chamber (MLR-351H, Sanyo, Osaka, Japan). To induce the salt and drought stress, respective solutions were added directly to the pots whereas foliar treatments were used for elicitor treatments (jasmonic acid and gibberellic acid). The samples were collected and stored at −80 • C until further use. Additionally, samples were also collected of plant organs (young leaves, mature leaves, and stems) and in vitro plant organs. For in vitro plant organs, nodal segments were surface sterilized using a previously established protocol [44] and inoculated on MS basal media (without any phytohormones). After 30 days, in vitro plant organs (young leaves, mature leaves, and stems) were collected and stored at −80 • C until further use.

Total RNA Extraction, cDNA Synthesis, and qRT-PCR Experiment Conditions
Total RNA was extracted from the plant tissues (±60-80 mg per sample) using the RNeasy Mini Kit (Qiagen, Hilden, Germany). The purity and integrity of the RNA samples were assessed using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA) and running the samples on a 1.2% agarose gel electrophoresis. For cDNA synthesis, TURBO DNA-free™ (Invitrogen, Waltham, MA, USA) Kit was used according to the manufacturer's instructions where 1µg (per sample) of the quality-checked gDNA-free RNA was used as the template. Ten commonly utilized reference genes in plant gene expression studies were selected and degenerated primers were designed (ATP-synthase, F1-ATPase, 60S rRNA, 28S rRNA, HLH71, β-tubulin, GST, GAPDH, EF-1α, and Ubiquitin). These genes were selected based on previous literature studies and the primers were designed based on the homologous sequence of related or other plant species using Primer3 software [(version 0.4.0) https://bioinfo.ut.ee/primer3-0.4.0/ (accessed on 15 August 2022)]. All primers designed were tested using a general PCR (C1000 thermocycler, Bio-Rad, Hercules, CA, USA) and confirmed on a 1.8% agarose gel electrophoresis. For the qRT-PCR analysis, CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories, Hercules, CA, USA) was used. The used reaction mixture was as follows: 5 µL of SYBR Green Master Mix (Applied Biosystems™, USA), 1 µL of primer mix (10 µM each), and 4 µL of cDNA (10 ng in total). The real-time thermocycler was programmed at an initial denaturation step at 95 • C for 10 min, followed by 40 cycles of 15 s at 95 • C and 1 min at 58.4 to 61 • C (based on the annealing temperature of the primer pairs). For the melting curves, stepwise heating was performed from 60 to 95 • C. To assess the reliability of the selected reference gene/s, the expression levels of phenylalanine ammonia lyase (PAL) were normalized using the best and the least stable gene under two elicitor treatments. The experiment was carried out with five biological and technical replicates.

Reference Gene Analysis
Gene expression stabilities of the candidate genes were examined by NormFinder [24], RefFinder (https://heartcure.com.au) (accessed on 15 August 2022) [25], comparative ∆Ct [22], BestKeeper [21], and geNorm [23] according to a previous study [12]. The NormFinder software relies on the stability values and the variation among them to rank the genes, and genes with the least variations are ranked as the best. The geNorm software package utilizes the expression stability value (M) to rank the candidate genes where the M value is inversely proportional to the gene ranking. BestKeeper relies on the standard deviation values of crossing point values (CP) or cycle threshold (Ct) and the coefficient of correlation (r) values. The RefFinder software integrates the ranking from other software to offer a complete ranking based on the geometric means of the ranking value.

Statistical Analysis
The relative expression of the PAL gene between the untreated control and the treated samples was compared using a two-sample t-test in Microsoft Excel 2021 with XLSTAT (a Statistical Software for Excel) (version 2022.1) (https://www.xlstat.com/en/) (accessed on 27 August 2022). The average Ct values from five biological replicates were used for the calculation of the relative PAL gene expression level. The gene expression levels were calculated using the 2 −∆∆Ct method [22].

Conclusions
In summary, the current study is the first systematic study to identify and validate the most suitable reference genes for qRT-PCR studies in M. officinalis or lemon balm. We ranked ten commonly used reference genes based on their expression stabilities. Based on the results, we recommend using a combination of EF-1α and GAPDH for the normalization of qRT-PCR in lemon balm. HLH71 was identified as the least stable gene. The identified optimal and the least stable genes were then validated with the PAL gene under elicitor stress. In conclusion, this study will provide a useful resource for more accurate and widespread experiments with qRT-PCR in M. officinalis and other related species.