Selection of Reference Genes for RT-qPCR Analysis in Coccinella septempunctata to Assess Un-intended Effects of RNAi Transgenic Plants

The development of genetically engineered plants that employ RNA interference (RNAi) to suppress invertebrate pests opens up new avenues for insect control. While this biotechnology shows tremendous promise, the potential for both non-target and off-target impacts, which likely manifest via altered mRNA expression in the exposed organisms, remains a major concern. One powerful tool for the analysis of these un-intended effects is reverse transcriptase-quantitative polymerase chain reaction, a technique for quantifying gene expression using a suite of reference genes for normalization. The seven-spotted ladybeetle Coccinella septempunctata, a commonly used predator in both classical and augmentative biological controls, is a model surrogate species used in the environmental risk assessment (ERA) of plant incorporated protectants (PIPs). Here, we assessed the suitability of eight reference gene candidates for the normalization and analysis of C. septempunctata v-ATPase A gene expression under both biotic and abiotic conditions. Five computational tools with distinct algorisms, geNorm, Normfinder, BestKeeper, the ΔCt method, and RefFinder, were used to evaluate the stability of these candidates. As a result, unique sets of reference genes were recommended, respectively, for experiments involving different developmental stages, tissues, and ingested dsRNAs. By providing a foundation for standardized RT-qPCR analysis in C. septempunctata, our work improves the accuracy and replicability of the ERA of PIPs involving RNAi transgenic plants.


INTRODUCTION
RNA interference (RNAi)-based genetically modified (GM) plants targeting insects have been developed and offer a new approach for insect control (Baum et al., 2007;Mao et al., 2007;Zha et al., 2011). Consuming transgenic maize that expresses an internal housekeeping gene, vacuolar ATPase subunit A, for example, significantly increases Diabrotica virgifera virgifera larval mortality (Baum et al., 2007). Technical and regulatory hurdles notwithstanding (Lundgren and Duan, 2013), commercialization of transgenic maize that expresses long dsRNAs for D. v. virgifera control appears likely in the near future (Kupferschmidt, 2013;Palli, 2014;Petrick et al., 2016;Tan et al., 2016). The impact of RNAi transgenic crops on non-target organisms (NTOs) is a major environmental concern. These NTOs include biological control agents that play essential role in integrated pest management (Romeis et al., 2008;Roberts et al., 2015;Xu et al., 2015).
Coccinella septempunctata (Coleoptera: Coccinellidae), the seven-spotted ladybeetle, is a generalist predator used for classical and augmentative biological control in many cropping systems. Its larvae and adults are voracious arthropod predators that also feed on pollen, nectar, and petals (Kalushkov and Hodek, 2004). C. septempunctata has been widely used as a 'model' NTO to evaluate the potential risks of Bacillus thuringiensis (Bt) transgenic crops (Harwood et al., 2005(Harwood et al., , 2007Alvarez-Alfageme et al., 2012). The mode of action of RNAi transgenic plants suggests that un-intended effects will likely occur via altered gene expression in NTOs (Berezikov, 2011), and reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) provides an important tool for detecting such changes. The importance of systematic criteria for selecting and evaluating reference genes used in RT-qPCR studies (Kalushkov and Hodek, 2004;Bustin et al., 2013;Li et al., 2013;Sinha and Smith, 2014;Zhu et al., 2014;Pan et al., 2015a;Wang et al., 2015) is illustrated by the fact that using one or multiple unsuitable reference genes for normalization can produce up to 20-fold differences in expression values and misinterpret gene expression results (Bustin et al., 2013;Hellemans and Vandesompele, 2014). Since RNAi-based insecticides and/or RNAi transgenic plants may cause lethal or sublethal effects in NTOs through changes in gene expression, RT-qPCR analysis without rigorous reference gene selection and validation may lead to inaccurate assessment of risks.
Reverse transcriptase-quantitative polymerase chain reaction is a powerful method for quantifying gene expression (Vandesompele et al., 2002). Although RT-qPCR is extensively used for measuring transcript abundance, variation in RNA extraction, RNA integrity and quality, enzymatic efficiency, and PCR efficiency can influence C q -values (Bustin et al., 2005;Strube et al., 2008). RT-qPCR generally involves normalization to the expression of a suite of appropriated reference genes in parallel. Even though reference gene transcript levels should ideally be stable across a range of different conditions, many commonly used reference genes differ dramatically across treatments (Kalushkov and Hodek, 2004;Bustin et al., 2013;Li et al., 2013;Sinha and Smith, 2014;Zhu et al., 2014;Pan et al., 2015a;Wang et al., 2015). Such high level of variation emphasizes the need to determine stable reference genes for RT-qPCR analyses on a case-by-case basis, even for the same species.
The goal of the current study is to select suitable reference genes for RT-qPCR analysis in C. septempunctata, specifically, for the normalization and analysis of vacuolar-type H + -ATPase (V-ATPase), a potential molecular target for RNAi transgenic plants.

Insects
Coccinella septempunctata (Coleoptera: Coccinellidae) pupae were collected from alfalfa at the north farm of University of Kentucky in June, 2015. Larvae and adults were reared in the laboratory at 23 ± 0.5 • C temperature, 16L: 8D photoperiod, and 50% relative humidity. They were provisioned with pea aphids (Acyrthosiphon pisum) reared in a greenhouse at 20-28 • C on fava bean, Vicia faba (Fabales, Fabaceae).

Biotic Factors
All developmental stages of C. septempunctata were sampled: eggs, all four larval instars (collected at the first day of each instar), pupae, female and male adults. Different body tissues, including the head, gut, and carcass (body without head or viscera), were dissected from larvae.

Abiotic Factor
For the dietary RNAi treatment, first-instar C. septempunctata larvae were supplied with 15% sugar solution containing one of the following: (1) in vitro synthesized dsRNAs from a target gene, C. septempunctata V-ATPase subunit A (dsCS) (dsCS forward: TA ATACGACTCACTATAGGGAGATCTCTTTTCCCATGT; dsCS reverse: TAATACGACTCACTATAGGGAGAGCATCTCGGCC AGAC); (2) a specific positive control, V-ATPase A 400 bp fragment amplified from D. v. virgifera; (3) a control gene, β-glucuronidase (dsGUS); and (4) a blank control, H 2 O (Yang et al., 2015). Neonate larvae emerging from their eggs were kept in separate petri dishes for 2 days; each larvae was provisioned on a daily basis with a 2 µl droplet containing 4 µg/µl dsRNA or the water control. On days 3, 5 individuals per treatment were collected for RT-qPCR analysis.
The number of sampled individuals per replicate in the developmental stage study was as follows: Egg stage: 15 eggs; first instar: five individuals; second instar: five individuals; third instar: three individuals; fourth instar: one individual; pupal: one pupa; adult male or female stage: one male or female individual. For the other biotic and abiotic factors, five individuals were sampled per replicate, and each experiment was replicated three times. Samples were placed in 1.5 ml centrifuge tubes, quickly frozen in liquid nitrogen, and stored at −80 • C prior to total RNA isolation.

Double-Stranded RNA Preparation
First-strand cDNA was prepared using 2.0 µg of total RNA with the M-MLV reverse transcription kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's recommendations. Pairwise comparison showed that the entire coding sequence of v-ATPase A from D. v. virgifera and C. septempunctata share a 81.0% nucleotide sequence similarity (Supplementary Figure  S1). The 400 bp region with the highest sequence similarity (85%) was selected as the template to synthesize arthropodactive dsRNAs (Supplementary Figure S2). A non-specific negative control, the β-glucuronidase (GUS) gene was cloned into pBTA2 vector and PCR amplified using gene specific primers, which amplified a 560 bp fragment containing T7 polymerase promoter region at the 5 end. PCR amplifications were performed in 50 µl reactions containing 10 µl 5 × PCR Buffer (Mg2+ Plus), 1.0 µl dNTP mix (10 mM of each nucleotide), 5.0 µl of each primer (10 µM each), and 0.25 µl of GoTaq (5 u/µl) (Promega). The PCR parameters were as follows: one cycle of 94 • C for 3 min; 35 cycles of 94 • C for 30 s, 59 • C for 45 s and 72 • C for 1 min; a final cycle of 72 • C for 10 min. The PCR product was used as template to generate dsRNA with the T7 MEGAscript kit (Ambion, Austin, TX, USA) following the manufacturer's protocol. The synthesized dsRNAs were suspended in nuclease-free H 2 O and quantified with a NanoDrop 2000c spectrophotometer and then stored at −20 • C.

Reference Gene Validation
For the tissue and dietary RNAi experiments, reference gene reliability was assessed by normalizing V-ATPase expression profiles with the two most-and two least-stable non-rRNA genes. For the development experiment, reference gene reliability was assessed by normalizing V-ATPase expression profiles with the three most-and three least-stable non-rRNA genes. Relative gene expression of V-ATPase was calculated using the 2 − Ct method (Livak and Schmittgen, 2001). One-way ANOVA was used to compare V-ATPase expression under each dietary RNAi treatments, across different developmental stages, and in different tissue types.

Candidate Gene Cloning and Performance
Four reference genes (Actin, ArgK, EF1A, and Tubulin) and one target gene (V-ATPase) were cloned based on degenerate primers. All candidate genes were expressed in C. septempunctata and visualized by a single amplicon of the expected size (Supplementary Figure S3). Gene-specific amplification of all candidate genes was confirmed by a single peak in meltcurve analysis (Supplementary Figure S4). The PCR efficiency (E), correlation coefficient (R 2 ), and linear regression equation characterizing each standard curve were given in Table 1. The standard curve of each gene was also provided (Supplementary Figure S5). The C q -values of all candidate genes under the three experimental conditions ranged from 10 to 29. The three ribosomal genes (18S, 28S, and 16S) showed mean C q -values less than 15 cycles. Actin, EF1A, Tubulin, NADH, and V-ATPase had C q -values ranging from 17 to 24 cycles. 28S and ArgK were the most-and least-expressed reference genes, respectively (Figure 1).

Expression Stability of the Reference Genes under Different Experimental Conditions
geNorm calculates the expression stability value 'M' for each reference gene. In the developmental stage study, 16S was the most stable gene. In different tissues, 18S and 28S were the most stable genes. In the dsRNA treatment, 18S and Tubulin were the most stable genes. Table 2 provides geNorm-based reference gene stability sequences for each experimental condition.
NormFinder calculates the expression stability value 'SV' for each reference gene. In the developmental stage study, 18S was the most stable gene. In different tissues, 28S and 18S were the most stable genes. In the dsRNA treatment, EF1A was the most stable gene. Table 2 provides NormFinder-based reference gene stability sequences for each experimental condition.
BestKeeper calculates the expression stability value 'SD' for each reference gene. For the developmental stage and tissue experiments, 16S was the most stable gene. For the dsRNA treatment, Actin was the most stable gene. Table 2 provides BestKeeper-based reference gene stability sequences for each experimental condition.
The C t method identifies potential reference genes by comparing expression of reference gene pairs within each sample. For the developmental stage and tissue experiments, 16S was the most stable gene. For the dsRNA treatment, Actin ranked as the most stable gene. Table 2 provides C t -based reference gene stability sequences for each experimental condition.

Optimal Number of Candidate Reference
Genes According to geNorm V-values across the different developmental stages, although never lower than 0.15, were lowest at V3/4. This implies that three reference genes were sufficient for normalization throughout developmental stages (Figure 3). The first V-value less than 0.15 emerged at V2/3 in both the tissue and dsRNA experiments, suggesting that two reference genes were sufficient for normalization (Figure 3).

Relative Gene Expression of V-ATPase
Among dsRNA treatments, V-ATPase expression differed when normalized to the two most-and least-stable non-rRNA reference genes (Figure 4). Expression of V-ATPase was most suppressed on day 3 in the dsDVV and dsCS treatments relative to the dsGUS and H 2 O controls (Figure 4).
Among developmental stages, V-ATPase expression differed when normalized to the three most-and least-stable non-rRNA reference genes. The expression level was different for each developmental stage under the two normalization conditions (Figure 5).
Among tissue types, V-ATPase expression was similar when normalized to the two most-and least-stable non-rRNA reference genes. When normalized to the most-and least stable reference genes, however, V-ATPase expression in the gut was about 7.7 and 22.4-fold higher, respectively, than in the carcass (Figure 6).

DISCUSSION
Although a wide range of techniques (e.g., cDNA microarray, subtractive hybridization, Western blot, Northern blot, RNA sequencing) can be used to study gene expression, the advantages of RT-qPCR in terms of its sensitivity and specificity, especially in non-model organisms, makes this an especially valuable tool. This technique, however, requires reference gene normalization to achieve reliable and comparable results (Kalushkov and Hodek, 2004;Li et al., 2013;Sinha and Smith, 2014;Zhu et al., 2014;Pan et al., 2015a;Wang et al., 2015). Our results confirm that the most stable reference genes can vary in different experimental conditions. While V-ATPase was the least stable C. septempunctata candidate gene in different tissues FIGURE 4 | Relative gene expression of the V-ATPase in C. septempunctata under dietary RNAi treatments. The relative gene expression levels of V-ATPase were normalized to the most suited non-rRNA (A, Actin and Tubulin) and the least suited (B, NADH and ArgK) 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 2 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.05). and dsRNA conditions, for example, it was very stable across different developmental stages. This is consistent with recent work suggesting that while V-ATPase was an unstable candidate gene for Hippodamia convergens (Coleoptera: Coccinellidae) in different tissue and dsRNA conditions, it was stable across developmental stages (Pan et al., 2015b).
The best-practice "Minimum Information for Publication of Quantitative Real-Time PCR Experiments" suggests using multiple reference genes in order to avoid biased normalization (Bustin et al., 2009). While two reference genes were adequate to analyze gene expression in different tissue and dietary RNAi conditions, we found that three reference genes were necessary across different developmental stages. The former result is consistent with previous work showing that two reference genes are sufficient for dependable gene normalization in different tissue and dietary RNAi conditions (Yang et al., 2015). The latter result, that more reference genes were required for FIGURE 5 | Relative gene expression of the V-ATPase in different developmental stages of C. septempunctata. The relative gene expression levels of V-ATPase in eggs (EG), first instar larvae (L1), second instar larvae (L2), third instar larvae (L3), fourth instar larvae (L4), pupae (P), male adults (MA), and female adults (FA) were normalized to the most suited non-rRNA (A, NADH, EF1A, and Tubulin) and the least suited (B, EF1A, Actin, and ArgK) reference genes, respectively. The relative expression level (fold) was calculated based on the value of the egg stage expression detected, which was assigned an arbitrary value of 1. Values are means ± SE. Different letters indicate significant expression differences among different developmental stages of C. septempunctata (P < 0.05).
developmental-stage analysis, likely reflects the fact that dramatic changes in gene expression occur during metamorphosis from one developmental stage to another (egg/larva, larva/pupa, pupa/adult). The C t -value of Actin is ∼24 at the egg stage, for example, and 16-19 at other developmental stages.
Levels of V-ATPase expression varied substantially following normalization to the most-and least-stable reference genes (Figures 4-6). This finding is consistent with previous work documenting condition-dependent variation in reference gene expression, and underlines the need for reference genes to be validated under particular experimental condition prior to their use (Bustin et al., 2013;Hellemans and Vandesompele, 2014;Yang et al., 2015).
Since rRNA makes up a large proportion of the total RNA pool (more than 80%), it is reflect by the three ribosomal genes (18S, 28S, and 16S) showed mean C q -values less than 15 cycles, FIGURE 6 | Relative gene expression of the V-ATPase in different tissues of C. septempunctata. The relative gene expression levels of V-ATPase were normalized to the most suited non-rRNA (A, NADH and Tubulin) and the least suited (B, Actin, and ArgK) reference genes, respectively. The relative expression level (fold) was calculated based on the value of the carcass expression detected, which was assigned an arbitrary value of 1. Values are means ± SE. Different letters indicate significant expression differences among different tissues of C. septempunctata (P < 0.05).
whereas mRNA makes up only 3-5%, thus, the use of rRNA for normalization of RT-qPCR may be problematic. Given this, it would have been preferable to include several mRNA species of the ribosomal machinery, e.g., 40S ribosomal protein S24, 40S ribosomal protein S18, 60S ribosomal protein L4 as opposed to only 18S, 28S, and 16S. Because many other studies still choose rRNA as the reference genes, however, we present the results of the three rRNA genes to illustrate how rRNAs are expressed under different experimental conditions in C. septempunctata. In addition, the least-expressed reference gene, like ArgK in C. septempunctata and C. maculata (Yang et al., 2015), is also not fit for the normalization.
This study, combined with our previous results, show that the three predatory ladybeetles, which share the same receiving environment (the maize field) and phylogenetically closely related to D. v. virgifera, are susceptive to ingested dsRNAs (Yang et al., 2015;Pan et al., 2015b). Additionally, dietary RNAi of dsDVV had no impacts on the other three NTOs including Apis mellifera, Sinella curviseta, and Danaus plexippus at both transcriptional and phenotypic levels (Pan et al., 2015a(Pan et al., , 2016Vélez et al., 2016). Our results are thus consistent with previous studies suggesting that the spectrum of dsRNA activity is expected to be narrow and species taxonomically related to the target organism are more likely to be susceptible (Baum et al., 2007;Bachman et al., 2013).
RNA interference has a wide range of applications in agriculture, especially for crop protection. Because many NTOs provide diverse ecosystem services (e.g., biological control, pollination, and decomposition), however, the effect of RNAi transgenic plants on species and the services they provide should be evaluated prior to commercialization. Our study thus provides a starting point for future work assessing the risks associated with RNAi transgenic plants. In addition, among different developmental stages and tissue types, V-ATPase expression differed when normalized to the two most-and least-stable non-rRNA reference genes, respectively. Expression of V-ATPase in the gut, for instance, ranged from 7.7 and 22.4-fold higher than in the carcass when normalized to the most-and least-stable sets of reference genes, respectively (Figure 6). Better accuracy in gene expression analysis not only can facilitate our investigation of gene function and evaluation of its efficacy for pest control, but also can improve the assessment of risks associated with RNAi transgenic plants.
Our results describe the selection and evaluation of stable C. septempunctata reference genes for use as internal controls in gene expression analysis. Eight endogenous reference genes were selected in 45 different samples for one abiotic (dietary RNAi) and two biotic (developmental stage and tissue type) conditions using five commonly used analytical methods. Different non-rRNA reference genes are recommended for each experimental condition: NADH, EF1A, and Tubulin across different development stages, NADH and Tubulin in different tissues, and Tubulin and Actin for the dietary RNAi experiment. Although the selection of reference genes is trait, species, condition and treatment specific, this study represents the first step toward establishing standardized RT-qPCR analysis in C. septempunctata. In addition, the methodology described here represents a critical step toward the development of an in vivo dietary RNAi toxicity assay for assessing the risks associated with RNAi transgenic plants.

AUTHOR CONTRIBUTIONS
XZ and HP conceived and designed research. HP, CY, and HZ conducted experiments. YL, LD, and XZ contributed reagents and analytical tools. HP and CY analyzed data. HP, EP, CY, and XZ wrote the manuscript.

FUNDING
The granting agencies have no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

ACKNOWLEDGMENTS
The authors are grateful to Hui Li for his assistance with the data analysis. This research was supported by a grant from USDA BRAG grant (3048108827), the National Natural Science Foundation of China (31501642), and a Special Fund for Agroscience Research in the Public Interest (201303028). The information reported in this paper (No. 16-08-026) is part of a project of the Kentucky Agricultural Experiment Station and is published with the approval of the Director.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.01672/ full#supplementary-material FIGURE 1 | Alignment of v-ATPase A ORFs between Coccinella septempunctata and Diabrotica virgifera virgifera. Identical nucleotides are highlighted in black boxes.
FIGURE 2 | The alignment of a highly conserved region within the ORFs of v-ATPase A from C. septempunctata and Diabrotica virgifera virgifera. This 400 bp fragment was selected as the target template to synthesis insecticidal dsRNAs.
FIGURE 4 | Melting curves of the eight candidate reference genes and one target gene in C. septempunctata.
FIGURE 5 | Standard curves of the eight candidate reference genes and one target gene in C. septempunctata.