Selection of reference genes for normalization of RT-qPCR data in gene expression studies in Anthonomus eugenii Cano (Coleoptera: Curculionidae)

The pepper weevil, Anthonomus eugenii Cano (Coleoptera: Curculionidae), is the main insect pest of peppers (Capsicum spp.) throughout the southern U.S. and a potential target for novel control methods that may require gene expression analyses. Careful selection of adequate reference genes to normalize RT-qPCR data is an important prerequisite for gene expression studies since the expression stability of reference genes can be affected by the experimental conditions leading to biased or erroneous results. The lack of studies on validation of reference genes for RT-qPCR analysis in A. eugenii limits the investigation of gene expression, therefore it is needed a systematic selection of suitable reference genes for data normalization. In the present study, three programs (BestKeeper, geNorm and NormFinder) were used to analyze the expression stability of candidate reference genes (β-ACT, ArgK, EF1-α, GAPDH, RPL12, RPS23, α-TUB, 18S and 28S) in A. eugenii under different experimental conditions. Our results revealed that the most stably expressed reference genes in A. eugenii varied according to the experimental condition evaluated: developmental stages (EF1-α, 18S and RPL12), sex (RPS23 and RPL12), low temperature (GAPDH and α-TUB), high temperature (α-TUB and RPS23), all temperatures (α-TUB and GAPDH), starvation (RPL12 and α-TUB), and dsRNA exposure (α-TUB and RPL12). Our study provides for the first time valuable information on appropriate reference genes that can be used in the analysis of gene expression by RT-qPCR in biological experiments involving A. eugenii.


Results
Primer specificity and efficiency of candidate reference genes. PCR products generated by each primer pair using cDNA from A. eugenii as a template were visualized as single bands of the expected size on 1.5% agarose gel ( Supplementary Fig. S1). The specificity of primer pairs was confirmed by sequencing of RT-PCR products and alignment with their corresponding gene fragment sequences. Additionally, the primer specificity was evaluated by melting curve analysis which showed the presence of a single peak ( Supplementary Fig. S2). A standard curve was generated for each primer pair using a serial dilution of the cDNA in order to calculate the correlation coefficient (R 2 ) and primer efficiency (E). E values varied from 93.65% to 108.39% and R 2 values were superior to 0.993 (Table 1). Expression profile of candidate reference genes. The mean quantification cycle (Cq) values varied considerably among the nine candidate reference genes, ranging from 8.82 (18S) in samples from different developmental stages to 28.83 (β-ACT) in samples from the starvation experiment (Fig. 1A,F). Overall, the candidate reference genes displayed similar expression patterns under different treatments. 18S and 28S had the lowest mean Cq values in all experimental conditions exhibiting the highest expression levels, whereas β-ACT and GAPDH showed the highest mean Cq values corresponding to the lowest expression levels ( Fig. 1A-G). The mean Cq values of the reference genes considering all treatments varied from 9.86 (18S) to 27.45 (GAPDH) and the standard deviations (SD) of Cq values ranged from 0.58 (RPS23) to 1.16 (β-ACT) (Fig. 1H).

Stability of candidate reference genes under different experimental conditions. Developmental
stages. EF1-α, 18S, 28S and α-TUB were determined to be the most stable reference genes by BestKeeper and NormFinder among different developmental stages, while RPL12, RPS23 and EF1-α were the top three most reliable reference genes according to geNorm. In contrast, geNorm and NormFinder ranked β-ACT and ArgK www.nature.com/scientificreports www.nature.com/scientificreports/ High temperature. 28S, 18S and α-TUB were ranked as the most suitable genes by BestKeeper, α-TUB, RPS23 and GAPDH by geNorm, and GAPDH, RPS23 and EF1-α by NormFinder under high temperature treatment. GeNorm and NormFinder identified 18S and 28S as the least stable genes, whereas BestKeeper indicated that β-ACT and ArgK had the worst performance ( Table 2).
All temperatures. When the gene stability was evaluated in insects submitted to high and low temperature stresses, GAPDH and α-TUB were the most stably expressed reference genes followed by ArgK or RPS23, according to geNorm and NormFinder. Based on BestKeeper, α-TUB, 28S and 18S displayed the most stable expression. β-ACT and EF1-α; 28S and 18S; β-ACT and 28S were the least stable genes according to BestKeeper, geNorm and NormFinder, respectively ( Table 2).
Starvation. BestKeeper analysis revealed that RPL12, RPS23 and α-TUB were the most stably expressed genes under starvation condition, while geNorm identified 18S, α-TUB and RPL12 as the most stable genes. According to results from NormFinder, RPL12, GAPDH and 18S showed the highest stability. All algorithms indicated that β-ACT and 28S were highly variable in their expression levels ( Table 2).
dsRNA. Based on geNorm and NormFinder, RPL12, RPS23 and α-TUB exhibited the highest expression stability for dsRNA exposure experiment, while BestKeeper indicated α-TUB, 28S and RPL12 as the most stable genes. BestKeeper and geNorm ranked GAPDH and β-ACT as the least stably expressed genes, while NormFinder indicated that 28S and β-ACT were the least stable (Table 2). overall ranking of candidate reference genes. As shown in Fig. 2 and Supplementary Table S1, the comprehensive ranking of candidate reference genes from the most to the least stable among the experimental conditions was as follows: optimal number of candidate reference genes. The pairwise variation (Vn/n + 1) analyses between two sequential normalization factors indicated that the pairwise variation V2/3 value was lower than the threshold value of 0.15 for sex, temperature, starvation and dsRNA treatments, suggesting that two reference genes are the optimal number of reference genes for accurate normalization of gene expression data in A. eugenii under these conditions. The pairwise variation V3/4 value was below the acceptable limit across developmental stages, therefore the use of three reference genes would be advisable to normalize the gene expression data (Fig. 3). The best combination of reference genes for normalization in each experimental condition is shown in Table 3.

Reference gene validation.
In order to validate some of the reference genes selected in this study, we evaluated the expression level of DNA-directed RNA polymerase II subunit RPB2 (RpII140) and Proteasome subunit alpha type2 (Prosα-2) genes in insects treated with the respective dsRNAs for these genes using the recommended set of reference genes (α-TUB and RPL12). We observed a significant gene knockdown of 92.2% and 96.5% in the insects injected with dsRNA targeting RpII140 and Prosα-2, respectively, compared to the control treatment in which the insects were injected with GFP dsRNA. However, only the insects fed on RpII140 dsRNA showed a significant decrease in gene expression level by 51.7% (Fig. 4).

Discussion
RT-qPCR has been used extensively to assess gene expression in entomological research [17][18][19][20] . Normalization of RT-qPCR data with reference genes is one of the most common strategies used to correct experimental errors introduced through the steps of RT-qPCR analysis; however, the choice of reference genes with low expression variation is necessary to guarantee valid normalization and avoid inaccurate gene expression quantification 21 www.nature.com/scientificreports www.nature.com/scientificreports/ Molecular studies provide important information about the genetic mechanisms underlying a variety of biological events and metabolic pathways, but only a few studies at the molecular level have been performed on A. eugenii. The selection of suitable reference genes for gene expression analysis will facilitate and boost such investigations in this important pest species.
In the present study, differences in the stability of potential reference genes were evaluated in order to select appropriate normalization factors for gene expression analysis in A. eugenii. Our results demonstrated that while some genes were ranked at the same position by BestKeeper, geNorm and NormFinder for a given condition, in general the stability ranking of the reference genes generated by these algorithms varied considerably. Variation in www.nature.com/scientificreports www.nature.com/scientificreports/ the ranking order has been observed in many studies and can be attributed to the different statistical approaches implemented in the algorithms 23,24 . To address this issue, a comprehensive ranking of reference genes was created based on the ranking value attributed by the algorithms geNorm, NormFinder and BestKeeper as performed previously [25][26][27][28] . Figure 3. Determination of the optimal number of reference genes for accurate normalization using pairwise variation (Vn/n + 1) analysis by geNorm.

Experimental condition Reference genes
Developmental stages  www.nature.com/scientificreports www.nature.com/scientificreports/ Our results indicated that α-TUB and RPL12 were consistently ranked as the most stable genes according to the overall ranking and at least one of these genes was included in the set of normalizer genes suggested for each experimental condition. α-TUB is the main component of microtubules that form the cytoskeleton structure, which plays an important role in several eukaryotic cellular processes, including cell mobility and division, and intracellular trafficking 29 . Our results are consistent with previous studies reporting that α-TUB was stable for developmental stages, sexes, tissues, temperature and photoperiod stresses in Empoasca onukii 8 , among various tissues from Coleomegilla maculata 30 and Mythimna separata 24 , for different temperatures in Phenacoccus solenopsis 31 as well as for developmental stages, tissues and sexes in Hermetia illucens 32 . The RPL12 gene encodes a structural protein of ribosomes and is involved in protein translation 33 . This gene exhibited high stability in Acyrthosiphon pisum exposed to temperature stresses 34 and in M. separata under different photoperiod and temperature treatments, and larval tissues 35 .
In contrast to our results, α-TUB was not a suitable reference gene for normalization across developmental stages of Cryptolestes ferrugineus 36 , in different tissues of Diaphorina citri 37 , developmental stages, sexes and in response to temperature stress in Propylea japonica 38 , and in nonviruliferous/viruliferous Frankliniella occidentalis 39 . These results indicate that the gene stability can be affected by the biotic and abiotic conditions and even the insect species evaluated. As a consequence, the selection of condition-specific reference genes is strongly recommended prior RT-qPCR analysis.
Our analyses further identified β-ACT and ArgK as the least stable reference genes for most of the experimental conditions indicating that these genes are not ideal for RT-qPCR data normalization in A. eugenii under the conditions tested. Although β-ACT and ArgK are traditional reference genes and have been used in many gene expression studies, they are not always stably expressed. Consistent with our results, many studies have shown unstable expression of β-ACT under variable conditions in a variety of species, including Bemisia tabaci 40 , M. separata 35 , F. occidentalis 39 , Henosepilachna vigintioctopunctata 41 , P. japonica 38 , Hippodamia convergens 42 and C. maculata 30 , as well as for ArgK in Spodoptera litura 43 , P. japonica 38 , C. maculata 30 and E. onukii 8 .
The evaluation of gene expression is a fundamental and routine analysis in RNAi-related studies, and therefore the choice of adequate reference genes is crucial to achieve precise results. We demonstrated that the selected reference genes were suitable to measure the relative expression of the target genes Prosα-2 and RpII140 in insects treated with dsRNA corresponding to these genes. The reference genes α-TUB and RPL12 identified in our work can be useful for normalization of target gene expression levels in further research on RNAi in A. eugenii. To date, neither the genomic or transcriptomic data are available for A. eugenii; however, this sequencing information could provide a valuable resource to gain a deeper understanding of different molecular mechanisms and to discover potential target genes that can be used in RNAi-mediated control methods against this insect pest. Future studies should focus on genome and transcriptome sequencing to fill the gap of genetic information on A. eugenii.
Despite the importance of accessing reference gene stability, these analyses are often not performed prior gene expression studies 44,45 . Our results highlighted the need to adopt this practice for proper normalization of RT-qPCR data. We found that the best combinations of reference genes that should be used as internal controls were EF1-α, 18S and RPL12 for developmental stages; RPS23 and RPL12 for sex, GAPDH and α-TUB for low temperature; α-TUB and RPS23 for high temperature; α-TUB and GAPDH for all temperatures; RPL12 and α-TUB for starvation; α-TUB and RPL12 for dsRNA treatment. The selected reference genes may be helpful in further gene expression studies in A. eugenii. In addition, the set of primers validated in this study can be used to evaluate the suitability of the candidate reference genes in experimental conditions other than those tested here.

Methods insect colony.
A. eugenii colony was maintained at 27 ± 10 °C, 30 ± 5% relative humidity, 14:10 h light:dark photoperiod and supplied with jalapeno peppers (Capsicum annuum) as an oviposition substrate and food source for both larvae and adults.
Reference gene selection, gene fragment cloning and primer design. Based on a literature search, we selected nine genes (β-ACT, ArgK, EF1-α, GAPDH, RPL12, RPS23, α-TUB, 18S and 28S) commonly used as reference genes and that have shown high stability in other insect species to investigate their suitability as reference genes for RT-qPCR in A. eugenii 7,8,30,34,43,46 . Degenerate primers for RPS23, 18S, 28S and RpII140 were manually designed based on conserved nucleotide sequence among other Coleoptera species. Primers for ArgK, EF1-α, RPL12, α-TUB and Prosα-2 were obtained from previous works 30,47 (Supplementary Table S2). PCR amplifications consisted of 5 µL of 10x PCR buffer, 8 µL of MgCl 2 (50 mM), 5 µL of dNTPs mix (10 mM of each nucleotide), 8 µL of each primer (10 µM), 1.25 µL of Taq DNA Polymerase (5U/µL) (ThermoFisher Scientific), 2 µL of cDNA (diluted 10×) and 12.75 µL of nuclease-free water. The PCR cycling conditions were as follows: one cycle of 95 °C for 3 min; 40 cycles of 95 °C for 30 sec, 50-55 °C for 30 sec and 72 °C for 30 sec; a final cycle of 72 °C for 5 min. PCR products were purified using QIAquick Gel Extraction Kit (Qiagen) and cloned into the pJET1.2/blunt cloning vector using the CloneJET PCR Cloning Kit (Thermo Scientific) following the manufacturer's protocol. Recombinant plasmids were transformed into One Shot ® TOP10 Chemically Competent Escherichia coli cells (Invitrogen) and sequenced by GENEWIZ company (South Plainfield, NJ, USA). According to the partial sequences of the genes, specific primers were designed using Primer3Plus software 48 (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi). Primers for β-ACT gene were designed based on the sequence available on NCBI (Accession number MH560343) and the primers for GAPDH were obtained from a previous work 49 (Table 1). experimental treatments. Developmental stages. Samples of each developmental stage of A. eugenii including eggs, first, second and third instars, pupae and adults were collected. Each biological replicate included 20 eggs, 18 first instar larvae, 12 second instar larvae, 10 third instar larvae, 6 pupae or 6 adults (3 females and 3 males).

Starvation.
A. eugenii adults were placed in a plastic vial and starved for 24 h at 27 ± 1 °C, 75 ± 5% relative humidity in a growth chamber. Each biological replicate consisted of six insects.
Temperature. To examine the effects of temperature on gene expression stability, A. eugenii adults were exposed to temperatures ranging from 5-40 °C at 5 °C increments for 3 h in a climate-controlled chamber. Insects maintained at 5-20 °C were included in the low temperature treatment, while insects maintained at 30-40 °C were included in the high temperature treatment. Insects exposed to 25 °C constituted the control group for both treatments. Each biological replicate consisted of six adults. Gene expression stability was also evaluated when insects exposed to all temperatures were taken into account.
dsRNA. For the RNAi experiments, we selected two target genes (RpII140 and Prosα-2) that could potentially be used in RNAi-mediated control methods against A. eugenii based on previous studies demonstrating mortality of insects exposed to dsRNA targeting homologs of these genes 50,51 . DNA template for the synthesis of dsRNA was amplified by PCR using gene-specific primers containing a T7 polymerase promoter sequence at the 5′ (Table 4) from plasmid DNA. The PCR reaction was purified using QIAquick Gel Extraction Kit (Qiagen) and the dsRNA synthesized and purified using the MEGAScript TM RNAi Kit (Invitrogen) according to the manufacturer's instructions. For the bioassay by microinjection, A. eugenii adults (11-12 insects) were injected dorsally with 0.5 µL of RpII140 dsRNA or Prosα-2 dsRNA at 1000 ng/µL using the IM-11-2 microinjector (Narishige). Control insects were injected with an equivalent amount of GFP dsRNA or nuclease-free water. After the injections, insects were maintained at 27 ± 1 °C, 75 ± 5% relative humidity and fed on pepper. Biological replicates of each treatment were collected three days post-injection. Each biological replicate consisted of one insect. In the feeding bioassay, a droplet consisting of 24 µL of 20% sucrose solution with green food dye containing RpII140 dsRNA or Prosα-2 dsRNA at a concentration of 500 ng/µL was offered to twelve A. eugenii adults placed in a plastic vial. GFP dsRNA and nuclease-free water were used as controls. The droplets were replaced on the third day. Biological replicates consisting of a pool of three insects were collected on the fifth day.
The collected samples from all experimental treatments performed in this study were placed in centrifuge tubes, rapidly flash-frozen in liquid nitrogen, and stored at -80 °C until RNA extraction. Each treatment included three biological replicates.
RNA extraction and cDNA synthesis. Total RNA was extracted using the RNeasy Mini Kit (Qiagen) and on-column genomic DNA digestion was performed using the RNase-free DNase Set (Qiagen) as recommended by the manufacturer. RNA samples were quantified by Nanodrop 1000 spectrophotometer (Thermo Scientific). The absence of DNA contamination and the RNA integrity was analyzed on a 1.5% agarose gel. Total RNA (500 ng) was reverse transcribed to cDNA using a High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems) following the manufacturer's protocol. The cDNA was diluted 50-fold with nuclease-free water for subsequent RT-qPCR assays.

Reverse-transcription quantitative pcR (Rt-qpcR).
The RT-qPCR assays were carried out using a BioRad CFX96 qPCR System (Bio-Rad) with an optical 96-well plate. Each RT-qPCR reaction mix contained 2 μL of cDNA diluted 50×, 5 μL of SsoAdvanced TM Universal SYBR ® Green Supermix (Bio-Rad), 0.2 μL of each primer at 10 μM and 3.6 μL of nuclease-free water for a total volume of 10 μL. Thermal cycling conditions were 95 °C for 2 min, followed by 40 cycles of 95 °C for 5 sec and 60 °C for 30 sec. To assure the specificity of the primers and to eliminate the possibility of primer dimer formation, melting curves ranging from 65 °C to 95 °C with 0.5 °C/5 sec increment were included after amplification. Non-template control (NTC) was used as a negative control for each master mix. Assays were performed with three biological replicates each comprising three technical replicates. The PCR amplification efficiency of the primer pairs was determined from the standard curve generated with 5-fold serial dilutions of cDNA.
Analysis of the stabilities of candidate reference genes. Three algorithms, geNorm accessed as part of the qbase+ analysis software from Biogazelle 52 (http://medgen.ugent.be/*jvdesomp/genorm/), NormFinder version 0.953 53 (https://www.moma.dk/normfinder-software) and BestKeeper version 1 54 (https://www. gene-quantification.de/bestkeeper.html) were used to evaluate the expression stability of the candidate reference genes. For NormFinder analyses, Cq values were transformed into non-normalized relative quantities according to the formula: (E) ΔCq where E represents the primer efficiency for each gene and ΔCq represents the lowest Cq value -Cq value of each sample. Raw Cq values were employed in geNorm and BestKeeper analyses. The comprehensive ranking of the reference genes was based on the geometric mean of geNorm, NormFinder, and BestKeeper results.

Primer
Sequence 5′−3′ Size (bp)  www.nature.com/scientificreports www.nature.com/scientificreports/ BestKeeper calculates the standard deviation (SD) and coefficient of variation (CV) based on the Cq values for each candidate reference gene. Genes with SD greater than 1 are considered unstable. It also estimates the BestKeeper Index and then calculates the correlation coefficient (r) between each candidate gene and the BestKeeper index. Genes with higher stability have higher r values. GeNorm evaluates the stability of the potential reference gene based on an "M" value which represents the average pairwise variation of a specific candidate reference gene with all other genes. The genes with the lowest M values are considered the most stably expressed. For homogeneous samples, suitable reference genes should have M values lower than 0.5, but if the samples are considered heterogeneous M values up to 1 are acceptable. GeNorm was used to determine the optimal number of reference genes required for accurate normalization. It calculates the pairwise variation values (Vn/n + 1) between two consecutively ranked normalization factors (NFn and NFn+1), where n is the number of reference genes used in the normalization factor. The stepwise inclusion of the subsequent more stable reference gene can result in an increase or decrease in Vn/n + 1 value and a Vn/n + 1 below 0.15 indicates that the inclusion of an additional reference gene is not necessary for normalization. NormFinder analyzes the intra-and inter-group variations to estimate the stability of each candidate reference gene. Lower stability values indicate higher gene expression stability.
Validation of the reference genes in A. eugenii. The combination of the most stable reference genes (α-TUB and RPL12) in insects treated with dsRNA was validated by evaluating the expression of Prosα-2 and RpII140 genes. The relative expression of the target genes was calculated using the 2 − ΔΔCt method 55 . Three technical and three biological replicates were performed in this analysis. Gene expression data were analyzed by one-way ANOVA followed by Tukey's HSD test at P < 0.05. All statistical analyses were performed using JMP Pro 13 Software (SAS Institute, Cary, NC). Data are presented as mean ± standard error of the mean (SE).

Data availability
All relevant data analyzed during this study are included in this article and its Supplementary Information files.