Optimization and identification of stable reference genes for RT-qPCR analysis in Capsicum annuum L.

Pepper (Capsicum annuum L.) is a globally significant vegetable crop, widely cultivated for use as both a culinary spice and fresh vegetable^[1]. The consumption of spicy peppers stimulates endorphin and dopamine release in the central nervous system, producing natural analgesic and euphoric effects, thereby making people feel happy^[2]. Beyond their culinary value, capsaicinoids demonstrate unique bioactive properties enabling multidisciplinary applications, including military, medical, maritime, and so on. However, the delayed initiation of whole genome sequencing in Capsicum species—combined with their complex genomic architecture (~3.5 Gb)—has significantly impeded functional genomics research. This directly restricts the molecular research and design breeding of pepper and the vigorous development of the pepper industry.

Since 2014, the progressive refinement of Capsicum genome assemblies has catalyzed a paradigm shift toward functional genomics research^[3−5]. In the study of the growth and development of peppers and their responses to biotic and abiotic stresses and under various hormone treatments, the expression analysis of related genes has become crucial. The level of gene expression often provides a reliable foundation for exploring gene functions^[6]. There are several conventional techniques for gene expression analysis, including real-time quantitative polymerase chain reaction (RT-qPCR), semiquantitative PCR, northern blotting, and in situ hybridization^[7]. Among them, RT-qPCR is renowned for its high degree of accuracy, robust specificity, and cost-effectiveness and has been widely used as a useful research tool for targeted gene expression analysis^[8−10]. The above benefits from the fact that when analysing the gene expression patterns detected by RT-qPCR technology, a stably expressed reference gene is often selected to correct and standardize the data, to eliminate the possible influence of factors such as RNA integrity, sample size, and reverse transcription efficiency on the RT-qPCR results^[11,12]. Furthermore, gene expression can also be affected by other factors such as tissue type, growth stage, and environmental conditions^[13]. This could have an impact on the precise investigation of gene function. Therefore, selecting appropriate reference genes with stable expression from thousands of genes plays a crucial role in ensuring the accuracy of RT-qPCR results for gene expression analysis^[14,15]. Common reference gene families, such as Ubiquitin, GAPDH, β-TUBULIN, Cyclophilin, β-actin, EF-1, etc., are frequently utilized as reference genes in numerous species^[16−18]. From 1900 to 2024, approximately 31,693 published pepper research articles were generated in the Web of Science database. Among them, about 2,693 articles were related to the analysis of pepper gene expression, and there were also some applications and reports of reference genes in RT-qPCR. In the published articles, research has demonstrated that the expression of UBI-3 and GADPH is the most stable in various tissues and that the stability of β-TUB and UBI-3 is strong under abiotic stress and hormone treatment. In contrast, UBI-1 and ACT were the most unstable and not suitable reference genes in the whole sample^[19]. In 2017, on the basis of the pepper transcriptome datasets of PGP (Pepper Genome Platform, http://passport.pepper.snu.ac.kr) and PGD (Pepper Genomics Database, http://ted.bti.cornell.edu/cgi-bin/pepper/index), CaREV05 and CaREV08 were used as the best reference genes for different development stages of pepper fruit^[20]. The stability of the published reference genes of peppers under different treatments varies to some extent. Therefore, to meet the various demands of pepper research, the stable reference genes still need to be further explored and optimized. In particular, based on the latest pepper transcriptome dataset and genomic information, the reference genes of pepper under different tissues, different hormone treatments, and different stresses have not yet been identified^[3−5].

In this study, based on the latest genome of pepper and the published datasets of pepper gene expression, potential reference genes were obtained. Their stability was comprehensively ranked by five different statistical algorithms. CaTubulin2, CaActin1, CaGAPDH1, and CaGAPDH2 were identified as relatively ideal reference genes. Further verification was conducted in hormone treatment, biotic and abiotic stress treatment, and published transcriptomes. It was found that CaActin1 and CaGAPDH2 are not only stable reference genes during pepper growth and development but also the optimal reference genes in hormone response and stress response. Overall, the above results provide genetic resources and a reference basis for future research on pepper gene expression.

The 'Zunla-1' peppers, whose genomic data have been reported, were selected as the experimental materials (provided by the pepper research group of the College of Horticulture, Hunan Agricultural University, Hunan, China). Three biological replicates per sample, each derived from a separate seedling, were analyzed, and this consistent methodology was applied to all subsequent experiments. Six parts, namely the roots, stems, true leaves of pepper at the six-leaf stage, the flowers and ovaries on the flowering day, and the fruits at the color-breaking stage, were selected to analyse the expression of reference genes in different organs. Another part of the material was taken from a greenhouse with a constant temperature of 22 °C, 12 h light/12 h dark, and a stable relative humidity of 65%–70%, and was used for hormone treatment and stress.

Hormone treatments and stress treatment

At the six-true-leaf stage of pepper seedlings, spray 100 μM GA₃^[21] or 100 μM IAA^[22], respectively. After 6 h of treatment, pepper leaves were collected and placed in an ultralow temperature refrigerator at −80 °C to analyse the stability of pepper reference genes under different hormone treatments.

At the six-true-leaf stage, waterlogging treatment was conducted, and root tissue samples were collected at 0 and 48 h after treatment. Pepper seedlings at the six-true-leaf stage were infected with Phytophthora capsici, and tissues were harvested 24 h after inoculation. Using materials from the same growth stage, plants were subjected to simulated natural drought conditions, with samples collected after 7 d of drought treatment.

RNA extraction and cDNA synthesis

The total RNA of the pepper plants was extracted with a Shanghai Promega Eastep Super total RNA extraction kit. Specifically, 30–50 mg pepper samples (leaves, roots, flowers, fruits, etc.) were used. The sample was ground with liquid ammonia in a mortar, placed in a sterile non-enzymatic centrifuge tube, and 300 μl of lysis solution, and vortexed to make pulp. Then, RNA was extracted from each sample according to the instructions of the RNA extraction kit of Promega.

A TIANGEN Fast King cDNA First Strand Synthesis Kit was used for synthesis. A 20 μl reaction system was established with 1 μg of total RNA, and the specific procedure was carried out according to the instructions.

Selection of candidate reference genes and primer design

Four commonly used reference families, namely, GAPDH, UIBIQUTIN, TUBULIN, and Actin, were selected on the basis of previously reported reference family genes. On the basis of the study of pepper reference genes in 2011^[19], combined with tomato genome data (SGN-VIGS [solgenomics.net]) and Arabidopsis gene data groups (TAIR-Home [Arabidopsis.org]), all genes associated with the GAPDH, UIBIQUTIN, TUBULIN, and Actin reference gene families and their related expression levels in various pepper tissues were obtained from the pepper website (Pepper Hub [hzau.edu.cn]). Reads per kilobase of transcript per million reads mapped (RPKM) value of each of the orthologous genes was collected. The average expression level and coefficient of variation of each gene in various pepper tissues were calculated. Screen out two to three candidate genes with an average expression level greater than 100 and an R² value less than 0.5 from each family. The CDSs of nine genes were downloaded from the pepper website (http://ted.bti.cornell.edu/cgi-bin/pepper/index Database (https://cornell.edu)). The primers used were designed on the fluorescence quantitative primer design website (Gen Script Pcr Primer Design) .

Detection of primer specificity and amplification efficiency of candidate reference genes

The cDNAs of six different tissue samples were mixed in equal amounts for primer-specific amplification detection. The PCR mixtures were as follows: cDNA (diluted two times), 1 μl each of the upstream and downstream primers (10 μM), 10 μL of Green Tag Mix (Nanjing Novartis), and 20 μL of ddH₂O. The PCR amplification procedure involved pre-annealing at 95 °C for 3 min; after denaturation at 95 °C for 15 s, annealing at 56 °C for 15 s, and extension at 72 °C for 1 min, 35 cycles of the above three steps were performed. The mixture was extended at 72 °C for 5 min and stored at 4 °C for later use.

The specificity of the PCR products was detected by 1.2% agarose gel electrophoresis after PCR. The steps of RT-qPCR were as follows: 2 μl of twofold diluted mixed cDNA was used as a template. According to the instructions of the ChamQ Universal SYBR qPCR Master Mix fluorescence quantitative kit, RT-qPCR was performed with a Roche LightCycler 96 fluorescence quantitative analyzer. After the reaction, melting curve analysis was performed to determine the specificity of the amplified product. The temperature was slowly increased from 60 to 95 °C. The fluorescence intensity of the sample was continuously measured to obtain the melting curve. If the melting curve was a single peak, the primer specificity was good.

After mixing the cDNA with an equal amount, it was diluted to zero, two, four, and six, respectively, as the templates for RT-qPCR. Three repeated experiments were conducted for each concentration. The Cq values of each candidate reference gene at different template dilutions were obtained via RT-qPCR. The standard curve was drawn with the log value of the template release multiple as the abscissa and the Cq value as the ordinate, and the slope (K) was calculated. The amplification efficiency (E) of the nine candidate reference genes was calculated via the formula E = 10⁻¹/K⁻¹.

Expression stability analysis of candidate reference genes in different tissues

In GeNorm, NormFinder, and BestKeeper software, the delta Ct value and average Cq value of each gene were used to analyse their expression stability. Finally, the stable and highly expressed candidate genes were screened via comprehensive evaluation.

Stability screening of candidate genes under different treatments

Further evaluate the top four stable candidate genes mentioned above. Taking the cDNA under different hormone treatments and different stress treatments as templates, further fluorescence quantitative expression analysis and determination were carried out to screen out the stable reference genes under different treatments.

Data analysis

The expression stability of the candidate reference genes was analysed via Genorm, NormFinder, and BestKeeper software. Genorm and NormFinder use the 2^−ΔΔCᴛ value for analysis; that is, the Cq value is converted to relative expression for analysis^[23]. Bestkeeper software is based on the input of the original data (Cq value) for analysis^[24]. Finally, the three software programs ranked the stability of each candidate gene according to the analysis results and selected the most suitable reference gene. The Cq values of the candidate genes under hormone and stress treatment in the later stage were analysed and mapped via GraphPad Prism 9.50 software.

In order to quickly understand the current status of pepper research, through the Web of Science database, the relevant literature on pepper and pepper gene expression analysis from 1988 to 2024 was retrieved. The corresponding number of documents for each year can be found in (Supplementary Table S1). Since 2007–2008, publications related to pepper have witnessed a rapid increase. Similar trends also exist in the publications related to gene expression among them (Fig. 1a). This might be the result of the synergy of the high-throughput sequencing technology iteration (the popularization of NGS (Next Generation Sequencing) represented by the Illumina platform), the integration of resources (International Program/Database), and the application pull around 2008, marking the entry of plant biology into a new era driven by big data^[25]. By reviewing the research in recent years, some commonly used reference genes under hormone treatment, biotic, and abiotic stress were revealed (Table 1). The Cq value of the reference gene should be around 20. Meanwhile, the difference in Cq values of the reference gene among samples should not exceed one to ensure the stability of the expression level of the reference gene. The stability of the above genes was verified through RT-qPCR. The results showed that only the UBI or Actin gene was relatively stable. More obviously, under disease and drought conditions, the Cq values of commonly used reference genes showed significant differences. Notably, the commonly used reference gene GAPDH exhibited significantly higher Cq values under hormonal treatment and various stress conditions, rendering it unsuitable as a reliable internal control under these experimental conditions (Fig. 1b). Therefore, it is urgently necessary to screen and optimize the reference genes in pepper, especially under stress conditions.

Figure 1. 

The existing status of reference genes in pepper. (a) Numbers of pepper gene expression publication from 1988 to 2024. (Based on the Web of Science database using the search word 'pepper' and 'pepper gene expression'). (b) Stability verification of published reference genes. (Significance analysis was performed using the two-tail Student's test (p < 0.05 indicated a significant difference).

Acquisition of nine candidate reference genes and specificity detection of their primers

According to previous studies, four evolutionarily conserved gene families (Tubulin, Actin, GAPDH, and Ubiquitin) were systematically selected for comparative analysis. The sequence information of these genes was obtained from relevant literature and databases for subsequent homology comparison analysis. A total of 87 genes (including nine GAPDH genes, 13 Ubiquitin genes, 27 TUBULIN genes and 38 Actin genes) were captured in the pepper genome website (Pepper Genomics Database, cornell.edu). The heatmap was plotted based on the FPKM values of these genes,along with their R² see (Supplementary Fig. S1) for details. Based on the FPKM values and coefficient of variation (CV) of each gene, genes with FPKM values ranging from 100 to 2,000, and a coefficient of determination (R²) ≤ 0.5 across all tissues, were selected. Only those genes exhibiting moderate expression levels in various tissues were considered candidate genes. A total of nine genes were selected and were respectively named CaGAPDH1, CaGAPDH2, CaUbiquitin1, CaUbiquitin2, CaTubulin1, CaTubulin2, CaActin1, CaActin2 and CaActin3 (Supplementary File 1), and design quantitative primers for each gene (Supplementary Table S2).

Specific primers were designed for amplification, and the specificity of the nine candidate reference genes was detected in 1.2% agarose gel electrophoresis using a mixture of cDNAs from each tissue as the template. The sizes of the obtained PCR products were all within the expected range and the bands were single (Fig. 2a), and the melting curves of the RT-qPCR primers showed a single peak (Fig. 2b). The amplification efficiency of each pair of primers was between 99% and 100%, and R² was greater than 0.99 (Fig. 2c). The results indicated that the primers designed for the above nine reference genes were all specific and effective, and could be used for further application and data analysis.

Figure 2. 

Specificity detection of nine pairs of primers. (a) The electrophoretic banding pattern of the amplification products of the nine reference genes. (b) The melting curve graphs of the nine reference genes. (c) Amplification efficiency of nine pairs of primers.

RNA extraction and quality assessment of different pepper tissues

Pepper seedlings at the six-leaf stage were selected as materials, and different tissues were used for RNA extraction (Fig. 3a). The extracted RNA was analysed by 1.2% agarose gel electrophoresis, which revealed clear bands with a 28 S to 18 S ratio of approximately 2:1, indicating no RNA degradation (Fig. 3b). The OD₂₆₀/OD₂₈₀ ratio of the obtained RNA samples ranged between 1.9 and 2.0 (Supplementary Table S3), indicating high RNA purity. These results demonstrate that the extracted RNA is suitable for subsequent experiments.

Figure 3. 

RNA samples and quality detection. (a) From left to right in turn is the six-leaf stage root, stems at six-leaf stage, six-leaf stage leaves, flowers on the day of anthesis, ovary at the flowering stage, fruits at the color-breaking stage. (b) The integrity of RNA was verified using 1.2% agarose gel electrophoresis, maintaining identical sample loading order to previous procedures.

Comparative stability analysis of candidate reference genes using different algorithms

The geNorm is a commonly used algorithm for determining gene expression stability using the M value. A lower M value indicates higher gene stability; an ideal reference gene should have an M value below 0.5 in homogeneous samples or below one in heterogeneous samples^[29]. The results revealed that the four genes with the greatest stability were CaTubulin2, CaActin1, CaGAPDH1, and CaGAPDH2 (Fig. 4a).

Figure 4. 

Expression stability value of nine candidate reference genes calculated by five complementary algorithms or evaluation parameters. (a) Genorm analysis of fluorescence quantitative data of nine genes. (b) Normfinder analysis of fluorescence quantitative data of nine genes. (c) Bestkeeper analysis of fluorescence quantitative data of nine genes. (d) ΔCt analysis of fluorescence quantitative data of nine genes. (e) The Cq values of fluorescence quantitative data of nine genes.

NormFinder calculated the intra-and inter-group expression variations of the genes and comprehensively evaluated their stability^[30]. The nine reference genes were ranked by their stability value (SV), where a lower SV corresponds to greater stability. The top four most stable genes identified were CaGAPDH2, CaUIBIQUTIN2, CaActin2, and CaTubulin2 (Fig. 4b). Additionally, BestKeeper analysis assesses reference gene stability by calculating Cq value variation (SD and CV), thereby selecting the most stably expressed genes least affected by experimental conditions^[24]. The CV and SD of the best reference gene should be lower. The top four genes were CaGAPDH1, CaGAPDH2, CaUIBIQUTIN1, and CaUIBIQUTIN2 (Fig. 4c). Furthermore, as a standard approach in qPCR analysis, the ΔCt method evaluates reference gene stability by measuring Cq value differences (ΔCt) between target and candidate reference genes. Decreasing SD values correlate with increasing stability of reference genes. The results of the ΔCT value analysis revealed that the top four genes related to the stability of the nine reference genes in different pepper tissues were CaUIBIQUTIN2, CaTubulin2, CaActin2, and CaGAPDH2 (Fig. 4d). However, the Cq values of each reference base were compared (Fig. 4e), and the reference genes with Cq values of 19–25 were selected^[31]. Integrating results from all five complementary algorithms or evaluation parameters, four candidate genes (CaTubulin2, CaActin1, CaGAPDH1, and CaGAPDH2) were prioritized for subsequent validation (Table 2).

Analysis of the stability of reference genes under GA₃ and IAA treatments.

Following the integrated ranking analysis, the four most stable candidates were evaluated under both GA₃ and IAA treatments. RT-qPCR quantification confirmed that CaTubulin2 and CaActin1 showed superior stability (Cq values < 25 cycles) under 100 μM GA₃ and IAA treatment, satisfying the stringent standards for reference gene selection (Fig. 5a, b).

Figure 5. 

Stability evaluation of four reference genes under hormonal and environmental stress. (a) The stability of four reference genes under IAA treatment. (b) The stability of four reference genes under GA3 treatment. (c) The stability of four reference genes under waterlogging treatment. (d) The stability of four reference genes under drought treatment. (e) The stability of four reference genes under disease treatment. Significance analysis was performedusing the two-tail Student's test (p < 0.05 indicated a significant difference).

Screening stable reference genes under biotic and abiotic stress treatments.

In parallel, the screening process was extended to assess the four candidate reference genes under waterlogging, drought, and disease conditions. RT-qPCR validation confirmed that CaGAPDH2 displayed exceptional expression stability (Cq < 25) during waterlogging and drought exposures (Fig. 5c, d), thereby qualifying as a reliable reference gene for these stress conditions. When pepper plants were infected with Phytophthora blight, both CaGAPDH1 and CaGAPDH2 exhibited excellent stability (Fig. 5e). However, CaGAPDH1 showed Cq values > 25, which exceeded the acceptable threshold for reference genes. Therefore, based on comprehensive evaluation, CaGAPDH2 was identified as the optimal reference gene for pepper studies under biotic and abiotic stress treatments.

Validation of expression stability for two reference genes in published pepper data

Expression profiles of CaActin1 and CaGAPDH2 across various tissues under different treatment conditions were obtained from the Pepper Hub website. Compared to untreated controls, CaActin1 expression levels remained highly stable at multiple time points following various hormone treatments (Fig. 6a). Similarly, CaGAPDH2 demonstrated remarkable stability under diverse stress conditions when compared to untreated plants (Fig. 6b). Meanwhile, analysis of the expression levels of the previously reported internal reference genes Actin and GAPDH in Pepper Hub (under different hormonal treatments in roots) revealed that their stability was relatively weaker than that observed in this study (Supplementary Fig. S2). Collectively, these data confirm that the reference genes screened in this study exhibit exceptional suitability—CaActin1 and CaGAPDH2 for normalizing gene expression data in pepper under diverse experimental conditions.

Figure 6. 

Analysis of expression stability of CaActin1 and CaGAPDH2 using published pepper data. (a) The expression of CaActin1 in roots under different hormone treatments. (b) The expression of CaGAPDH2 in roots under different stresses.

Quantitative real-time PCR (qRT-PCR) is the most widely used technique for gene expression analysis, with reference genes serving as critical controls for data normalization^{[16, 32]}. Appropriate reference genes can effectively correct the errors caused by RNA quality, reverse transcription efficiency, and operation technology, which is conducive to obtaining real expression differences in target genes^[33,34]. Extensive cross-species studies have identified commonly used reference gene families such as Actin, Tubulin, Ubiquitin, GAPDH, and EF1-α^[35,36]. However, their expression profiles vary significantly across varieties, tissue types, and experimental conditions, necessitating context-specific selection and optimization of reference genes^[37,38].

As a globally traded horticultural commodity, pepper contributes substantially to the agricultural economy through its diversified uses in food, medicinal, and defense applications^[39,40]. To date, the reference genes of pepper under diverse treatment conditions remain incompletely characterized, and some proposed reference genes require re-evaluation based on the updated pepper genome assembly (Fig. 1).

In this study, nine candidate genes were screened from four classical gene families (Actin, Tubulin, Ubiquitin, and GAPDH) based on the newly assembled pepper genome expression data (Supplementary Fig. S1)^[41]. The nine genes were subsequently preliminarily screened via RT-qPCR using different tissues of Zunla-1 varieties as templates. Using the melting curve and Cq values as the basis, three software tools (GeNorm, Norm Finder, and BestKeeper) were used to analyse the experimental data. Discrepancies were observed in the identification of optimal reference genes when applying five complementary algorithms or evaluation parameters, likely attributable to the distinct computational methodologies and weighting strategies employed by each analytical approach^[42]. Using a variety of algorithms to evaluate the experimental results can reduce the one-sidedness of the experimental results and increase the accuracy of the experimental results^[19]. Through integrated analysis of results from three distinct algorithms combined with ΔCt and Cq values across all tissue types, four optimal candidate reference genes were identified: CaTubulin2, CaActin1, CaGAPDH1, and CaGAPDH2 (Fig. 4, Table 2). These genes were subsequently validated for their expression stability under various experimental treatments. Among the above four genes, CaTubulin2, presented the highest stability in different tissues, which was consistent with the results reported by Zhao et al.^[19]. Notably, CaGAPDH2 also exhibited considerable stability in multiple tissues in our current study (Fig. 5).

Phytohormones play pivotal roles throughout plant growth and development, being essential for maintaining physiological balance^[43]. Research in this field encompasses numerous genes involved in biosynthesis, metabolism, and signal transduction pathways^[44,45]. When investigating gene expression within these regulatory networks, reliable reference genes become particularly essential. Among these, auxin and gibberellin (GA) represent two phytohormones critical for pepper growth and development, while also mediating stress responses^[46−49]. Building upon the four tissue-stable reference genes identified previously, hormonal treatments were conducted with 100 μM IAA and 100 μM GA₃ on pepper tissues to validate reference gene stability under phytohormone treatment^[21]. By analysing the data of RT-qPCR results, CaActin1 showed high stability under different hormone treatments (Figs 5a, b, and 6a).

In addition, drought and waterlogging, as common abiotic stresses in plants, also have an important impact on the growth and development of pepper. Long-term flooding (LT-WL) seriously affects and limits the growth and distribution of pepper^[50,51]. Pepper is in the LT-WL state, its aerobic respiration is inhibited, and anaerobic respiration is increased, which not only weakens energy metabolism, but also produces harmful substances, which have a great impact on its production and economic yield.^[52,53]. In this study, building upon our initial screening of tissue-stable reference genes, CaGAPDH2 was further identified as the most reliable reference gene under both waterlogging and drought stress conditions from among the four candidate genes (Fig. 5c, d). Furthermore, biotic stress represents a critical research focus in plant science. Phytophthora capsici, a devastating soil-borne pathogen, poses severe threats to pepper cultivation. The oomycete pathogen can cause severe yield losses—reaching up to 100% under favorable conditions—leading to substantial economic impacts^[54,55]. Due to the lack of effective prevention and treatment measures, it is very important to explore the regulatory mechanism of immune genes of pepper blight^[56]. Therefore, the most stable reference gene among the four candidate genes under Phytophthora inoculation was further explored. Through experimental validation and data analysis, CaGAPDH2 was found to be the most stable and highly expressed gene (Figs 5e and 6b). In general, CaGAPDH2 was identified as a relatively stable reference gene under both biotic and abiotic stresses. However, Zhao Jing et al. screened β-TUB and UBI-3 as more stable reference genes under abiotic stress in 2011^[19]. The above results are different from those of previous studies, which may be due to the updated pepper genome, different varieties, different hormones, and different stress treatments.

Finally, a comprehensive evaluation of candidate reference genes was conducted using recently published expression databases (Fig. 6). Analysis of gene expression patterns under various hormonal treatments and stress conditions revealed that the selected genes (CaActin1 and CaGAPDH2) exhibited robust expression profiles across different treatment conditions and time points. These findings underscore the critical importance of continually updating and optimizing reference gene selection for accurate gene expression studies.

Under phytohormone treatments, CaActin1 demonstrated optimal stability as a reference gene, while CaGAPDH2 proved to be the most stable under both biotic and abiotic stress conditions in pepper.