Selection of reference genes for RT-qPCR analysis in a predatory biological control agent, Coleomegilla maculata (Coleoptera: Coccinellidae)

Reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) is a reliable technique for quantifying gene expression across various biological processes, of which requires a set of suited reference genes to normalize the expression data. Coleomegilla maculata (Coleoptera: Coccinellidae), is one of the most extensively used biological control agents in the field to manage arthropod pest species. In this study, expression profiles of 16 housekeeping genes selected from C. maculata were cloned and investigated. The performance of these candidates as endogenous controls under specific experimental conditions was evaluated by dedicated algorithms, including geNorm, Normfinder, BestKeeper, and ΔCt method. In addition, RefFinder, a comprehensive platform integrating all the above-mentioned algorithms, ranked the overall stability of these candidate genes. As a result, various sets of suitable reference genes were recommended specifically for experiments involving different tissues, developmental stages, sex, and C. maculate larvae treated with dietary double stranded RNA. This study represents the critical first step to establish a standardized RT-qPCR protocol for the functional genomics research in a ladybeetle C. maculate. Furthermore, it lays the foundation for conducting ecological risk assessment of RNAi-based gene silencing biotechnologies on non-target organisms; in this case, a key predatory biological control agent.


Results
Performance of RT-qPCR primers. All gene candidates tested were visualized as a single amplicon of expected size on a 2.0% agarose gel ( Figure S1). Furthermore, gene-specific amplification was confirmed by a single peak in the melting-curve analysis (Fig. 1). The linear regression equation, correlation coefficient, and PCR efficiency for each standard curve are shown in Table 1. Additionally, the standard curve of each gene is shown in Figure S2.
C t values of candidate reference genes. The C t values of these 16 candidate reference genes under the four experimental conditions ranged between 9 and 35. The average C t value of the four ribosomal genes, including 18S, 28S, 12S, and 16S, was under 15 cycles. Actin and NADH showed an averaged C t value of less than 20 cycles. The averaged C t values of EF1A, GAPDH, Tubulin, RPS24, HSP70, HSP90, RPS18, RPL4, and V-ATPase were between 20 and 25 cycles. 18S and ArgK were the most and the least expressed reference gene, respectively (Fig. 2).
Stability of candidate reference genes under specific experimental conditions. Developmental stages included eggs, all four larval instars (collected at the first day of each instar), pupae, adult females and males. Tissues, including head, gut, and carcass, were dissected from C. maculata larvae of various instars. For the sex, gene expression profiles were, respectively, investigated in adult females and males. For dietary RNAi study, four dietary treatments were included; artificial diets containing dsRNAs from dsDVV, dsCM, dsGUS, and H 2 O (vehicle control). The average expression stability value (M-value) is used by geNorm to determine the best set of reference genes. Recommended M values for geNorm are M < 0.5 for homogeneous samples and M < 1 for heterogenous samples. Here, the lower the M-value coefficient, the higher the stability ranking. Developmental stage analyses showed RPS24 and RPS18 were co-ranked as the most stable genes. Tissue-specific experiments indicated that Tubulin and GAPDH were the most stable genes. Sex results showed that HSP70 and RPS24 were co-ranked as the most stable genes. Dietary RNAi treatment revealed that 12S and 18S were the most stable genes. Table 2 shows the overall ranking of these reference gene candidates from the most-to-least stable ones under each experimental condition.
A low stability value (SV) suggests a more stable gene by NormFinder. For the developmental stage experiment, V-ATPase was the most stable gene. Tissue-specific experiments indicated that 12S was the most stable gene. Sex results showed that 16S was the most stable gene. The 18S gene was considered the most stable for the dietary RNAi treatment experiment. The overall order based on NormFinder from the most-to-least stable reference genes is shown in Table 2.
The stability of a gene is inversely proportional to the standard deviation (SD) value as computed by BestKeeper program. Those with SD > 1 are excluded. EF1A was determined to be the most stable gene for the developmental stage experiment, compared to the tissue experiment where 16S was considered to be the most stable. GAPDH was the most stable gene for both sexes. 18S was shown to be the most stable gene for RNAi experiments. The overall order based on BestKeeper from the most-to-least stable reference genes are also found in Table 2.
The Δ C t method depends on a concept similar to that of geNorm, it also relies on relative pair-wise comparisons. Using raw C t values, the average SD of each gene set is inversely proportional to its stability. Here, V-ATPase was the most stable gene for the developmental stage experiment, while to the tissue-specific experiments, where was shown 18S to be the most stable gene. 16S was the most stable gene for both sexes and 18S was the most stable gene for RNAi experiments. The overall order based on the Δ C t method, from the most-to-least stable reference genes is shown in Table 2.
Comprehensive ranking of reference genes. RefFinder is a comprehensive program that integrates all four above-mentioned software tools to rank the candidate reference genes based on their stability. The following rankings are listed in order of most-to-least stable reference genes. For the developmental stages, the comprehensive ranking was V-ATPase, RPS18 Quantitative analysis of candidate reference genes based on geNorm. Each experimental condition may demand a different set of requirements for normalizing the RT-qPCR data. The first V-value < 0.15 emerged at V5/6, suggesting that five reference genes are needed for reliable normalization throughout developmental stages (Fig. 4). In regard to tissue-specific and dietary RNAi experiments, the first V-value < 0.15 emerged at V2/3, suggesting that two reference genes are necessary for the reliable normalization (Fig. 4). Based on the same principle, three reference genes are required for the reliable normalization of ladybeetle samples with different sex as the first V-value < 0.15 appeared at V3/4 (Fig. 4).  gene instead of the reference gene, which reflected by the highly varied expression levels under the dietary RNAi treatments (Fig. 3D).

Discussion
Housekeeping genes, constitutively expressed to maintain basic cellular functions, are the conventional choice for a standardized reference 33 . Interestingly, there is, in fact, no "universal" reference gene that is stably expressed and applicable for all cell and tissue types across various experimental conditions 42,[44][45][46][47][48][49] . Therefore, each candidate reference gene should be evaluated under specific experimental conditions 42,50 . Our results demonstrate that the suitable reference genes can be different in response to diverse biotic and abiotic conditions (Table 2; Fig. 3). For example, GAPDH was stably expressed in C. maculata under the tissue-and sex-specific conditions; however, its expression was highly variable among different developmental stages. This is consistent with the results from the convergens ladybeetle, Hippodamia convergens (Coleoptera: Coccinellidae), in which the expression of GAPDH was stable among different tissue types and sexes, but variable across different developmental stages 45 .
RT-qPCR is arguably the most widely used molecular technique for the detection and quantification of nucleic acids 50 . However, it is far from being a "gold standard" because of the lack of transparency, standardization and technical/quality controls 38 . Hellemans and Vandesompele 39 estimated the average difference in expression level of a gene of interest after normalization with any of two randomly selected non-validated reference genes is between 3 and 6-fold among 10-25% of the case studies. Such inconsistency makes it impossible to draw a conclusion with biological or clinical relevance. To avoid biased normalization, more and more researchers have started to embrace the idea of using multiple reference genes to analyze gene expression 42,[44][45][46][47][48][49] .
Determination of the optimal number of reference genes usually produces a trade-off between accuracy and practicality. In this study, five reference genes are required for reliable normalization under different developmental stages. In comparison, no more than three reference genes were required for reliable normalization under different sex, tissue types and dietary RNAi treatments. Metamorphosis has significant impact on the cellularity and consequently gene expression across the developmental stage. For examples, the C t value of ArgK was approximately 27 from egg to the fourth instar larva, whereas C t value increased to 35 at pupa and adult stage. Similarly, GAPDH had a C t value of 27 at the pupa stage, whereas it decreased to 23 at the other stages.
Our analyses demonstrate a dynamic shift in gene expression levels when normalized to reference genes that were determined to be the most and least suitable for a given treatment conditions (Fig. 5A,B). This provides a case-specific framework for selecting the most appropriate genes for normalization, as comparative measurements can yield varying results when using different gene sets to normalize data. Our study is consistent with previous studies showing how the variability in reference gene expression under variable experimental conditions can statistically affects study outcomes, thus strongly supporting the argument for reference gene validation prior to their use experimentally [51][52][53] .
The mRNA expression level of V-ATPase in C. maculata was apparently affected by dietary RNAi treatments. V-ATPase expression was significantly reduced under the dsDVV and dsCM treatments compared to the dsGUS and H 2 O controls (Fig. 5A). Coleomegilla maculate, a conventional NTO surrogate species which serves as a biological control agent, seems to be susceptive to a systemic exposure to the ingested dsRNAs. As a sequence-specific gene silencing tool, RNAi has a great potential in agricultural applications, either through crop improvements or pest/disease controls. Before this novel pest control strategy can be regulated/commercialized, the ecological risk assessment of RNAi-based controls on NTOs must be preceded. Our study provides a road map for future investigations on the risk assessment of RNAi-based gene silencing biotechnologies, including RNAi insecticides and transgenic RNAi crops.
In summary, expression profiles of 16 candidate reference genes under four experimental conditions (different tissue types, developmental stages, sex, and dietary RNAi) were investigated using five readily available algorithms (geNorm, NormFinder, BestKeeper, Δ C t method, and RefFinder). A suite of reference genes were specifically recommended for each experimental condition. These combined results reaffirm that there is no single universal reference gene suitable for all conditions, and reference genes can respond differently to various experimental conditions. This study represents the critical first step to establish a standardized RT-qPCR protocol for the functional genomics research in a ladybeetle C. maculate. Furthermore, it lays the foundation for conducting ecological risk assessment of RNAi-based gene silencing biotechnologies on non-target organisms; in this case, a key predatory biological control agent.

Materials and Methods
Insect cultures. Coleomegilla maculata (Coleoptera: Coccinellidae) was collected from cardoon, Cynara cardunculus, at the University of Kentucky in August, 2014. Larvae and adults were maintained in the laboratory and provisioned with pea aphids, Acyrthosiphon pisum, at 23 ± 0.5 °C, 16L: 8D photoperiod, and 50% relative humidity. Pea aphid clones were kindly provided by Dr. John Obrycki (University of Kentucky), and were maintained at 20-28 °C on fava bean seedlings, Vicia faba (Fabales, Fabaceae), in a greenhouse.  Experimental conditions. Biotic factor. The different developmental stages included eggs, all four larval instars (collected at the first day of each instar), pupae, and adults (including both females and males). Tissue types, including head, gut, and carcass (the remaining tissues that removed head and viscera) were dissected from various instars of C. maculate larvae. For different sex, one adult female and male were collected, respectively.
Abiotic factor. For dietary RNAi treatments, the first-instar larvae were fed with an artificial diet containing 15% sucrose solution mixed with chemically synthesized dsRNAs from 1) a target species, the western corn rootworm, D. v. virgifera (dsDVV, Forward: TAATACGACTCACTATAGGGAGAGCTCTTTTCCCATGTGTAC; Reverse: TAATACGACTCACTATAGGGAGAGCATTTCAGCCAAACG), and 2) a NTO, C. maculate (dsCM, Forward: TAATACGACTCACTATAGGGAGATCTCTTTTCCCATGT; Reverse: TAATACGACTCACTATAGGGAGAGCATCTCGGCCAGAC). The molecular target here is V-ATPase subunit A, an energy related housekeeping gene. Controls included an exogenous control gene β-glucuronidase from bacteria (dsGUS, Forward: TAATACGACTCACTATAGGGAGAGGGCGAACAGTTCCTGATTA; Reverse: TAATACGACTCACTATAGGGAGAGGCACAGCACATCAAAGAGA), and H 2 O, the vehicle control. At the beginning of the experiment, C. maculata neonates that hatched in less than 24 hours were kept individually in each petri dish. Each neonate was provisioned with a 2 μ l droplet containing 1 μ l of dsRNA (8 μ g/μ l) and 1 μ l of 30%     RPL4 and HSP70) reference genes, respectively. For dietary RNAi, ladybeetle larvae were exposed to an artificial diet containing 15% sugar solution and 4.0 μ g/μ l dsRNAs for two days (see Materials and Methods for details). The transcript levels of V-ATPase in newly emerged (0 day) untreated larvae were set to 1, and the relative mRNA expression levels in dsRNA-fed larvae were determined with respect to the controls. Values are means ± SE. Different letters indicate significant differences between the treatments and controls (P < 0.01). sucrose solution on a daily basis. For the first two days, a total of 16 μ g of dsRNA were provided to each neonate. On day-3, five individuals from each treatment were collected as one sample for the subsequent RT-qPCR analysis. For the developmental stage, a total of 15 eggs were collected as one biological replicate, while one pupa was collected, individually, as one replicate. For the remaining developmental stages, and all other biotic and abiotic conditions, approximately five individuals were collected for each treatment, and each experiment was repeated three times independently. All collected samples were flash frozen in liquid nitrogen and stored at − 80 °C in 1.5 ml centrifuge tubes. All the experiments were conducted at 23 °C with a photoperiod of 16: 8 (L: D).
Candidate reference genes and primer design. A total of 16 candidate reference genes commonly used in RT-qPCR analyses in other insect species were selected (Table 1). Primers for 12S, 16S, 18S, and 28S were designed based on the sequences obtained from NCBI. For the other seven genes including Tubulin, RPS24, HSP70, HSP90, NADH, RPS18, and RPL4 genes, primers were designed based on the sequences from a transcriptome of C. maculate 43 (Table S2). For the ArgK, EF1A, GAPDH, Actin, and V-ATPase genes, degenerate primers were designed using CODEHOP (http://blocks.fhcrc.org/codehop.html) according to conserved amino acid residues among Coleoptera species (Table S1). Conditions for PCR amplifications have been described previously 44,45 . PCR products were cloned into the pCR4-TOPO vector (Invitrogen, Carlsbad, CA), and sequenced. After the identities of these reference genes were confirmed (Table S2), primers for the subsequent RT-qPCR analyses were designed online, https://www.idtdna.com/Primerquest/Home/Index.

Reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR). The information regard-
ing RT-qPCR analysis has been described previously 44,45 . In brief, gene-specific primers (Table 1)  Data analysis. One way ANOVA was used to compare the gene expression of V-ATPase under each dietary RNAi treatments. Stability of the 16 candidate reference genes were evaluated by algorithms geNorm 33 , NormFinder 54 , BestKeeper 55 , and the Δ C t method 56 . Finally, RefFinder (http://www.leonxie.com/referencegene. php), a comprehensive software platform integrating all four algorithms, provided an overall ranking of the stability/suitability of these candidates 57 . Pairwise variation (V), as determined by geNorm, is an index for determining the optimal number of reference genes for accurate RT-qPCR normalization. A cut-off value for pairwise variation of 0.15 was recommended by Vandesompele et al. (2002) 33 . Beginning with two genes, this algorithm continuously adds another gene and recalculates the normalization factor ratio. If the added gene does not increase the normalization factor ratio over the proposed 0.15 cut-off value, the starting pair of genes is considered sufficient for normalizing data, otherwise, more genes should be incorporated.