Front. Plant Sci. Frontiers in Plant Science Front. Plant Sci. 1664-462X Frontiers Media S.A. 10.3389/fpls.2016.01672 Plant Science Original Research Selection of Reference Genes for RT-qPCR Analysis in Coccinella septempunctata to Assess Un-intended Effects of RNAi Transgenic Plants Yang Chunxiao 1 2 3 Preisser Evan L. 4 Zhang Hongjun 3 5 Liu Yong 2 Dai Liangying 1 Pan Huipeng 3 6 * Zhou Xuguo 3 * 1College of Plant Protection, Hunan Agricultural University Hunan, China 2Institute of Plant Protection, Hunan Academy of Agricultural Sciences Hunan, China 3Department of Entomology, University of Kentucky, Lexington KY, USA 4Department of Biological Sciences, University of Rhode Island, Kingston RI, USA 5Institute for the Control of Agrochemicals, Ministry of Agriculture Beijing, China 6Department of Entomology, South China Agricultural University, Key Laboratory of Bio-Pesticide Innovation and Application of Guangdong Province Guangzhou, China

Edited by: Joachim Hermann Schiemann, Julius Kühn-Institut, Germany

Reviewed by: Joe Hull, Agricultural Research Service (USDA), USA; Detlef Bartsch, Bundesamt für Verbraucherschutz und Lebensmittelsicherheit, Germany

*Correspondence: Xuguo Zhou, xuguozhou@uky.edu Huipeng Pan, hppan0623@gmail.com

This article was submitted to Plant Biotechnology, a section of the journal Frontiers in Plant Science

 08 11 2016 2016 7 1672 28 08 2016 24 10 2016 Copyright © 2016 Yang, Preisser, Zhang, Liu, Dai, Pan and Zhou. 2016 Yang, Preisser, Zhang, Liu, Dai, Pan and Zhou

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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.

 Coccinella septempunctata RT-qPCR reference gene RNAi transgenic plants environmental risk assessment plant incorporated protectant 3048108827 U.S. Department of Agriculture10.13039/100000199 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, 2007; Alvarez-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 Cq-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. A total of eight candidate reference genes were investigated: NADH dehydrogenase subunit 5 (NADH), elongation factor 1α (EF1A), β-actin (Actin), α-tubulin (Tubulin), arginine kinase (ArgK), 28S ribosomal RNA (28S), 16S ribosomal RNA (16S), and 18S ribosomal RNA (18S). The consistency of each reference gene was evaluated under one abiotic (dietary RNAi) and two biotic (developmental stage and tissue) conditions. Expression of the target gene, V-ATPase, was investigated under each of the three treatments using different normalization strategies.

 Materials and Methods 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: TAATACGACTCACTATAGGGAGATCTCTTTTCCCATGT; dsCS reverse: TAATACGACTCACTATAGGGAGAGCATCTCGGCCAGAC); (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, H2O (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.

 Total RNA Extraction and cDNA Synthesis

Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) in accordance with previously published methods (Yang et al., 2014, 2015). DNase treated total RNA was denatured at 75°C for 5 min and immediately chilled on ice. The concentration of RNA was determined using a NanoDrop 2000c Spectrophotometer. RNA concentrations were as follows: 287.9 ± 87.0 ng/μl [mean ± standard error of the mean (SEM)] for eggs, 392.0 ± 45.1 ng/μl for the first instars, 854.0 ± 62.2 ng/μl for the second instars, 480.4 ± 12.9 ng/μl for the third instars, 879.2 ± 153.8 ng/μl for the fourth instars, 727.8 ± 147.1 ng/μl for pupae, 622.6 ± 109.1 ng/μl for male adults, 558.1 ± 78.7 ng/μl for female adults, 501.8 ± 72.5 ng/μl for heads, 1095.5 ± 39.4 ng/μl for carcasses, 597.6 ± 62.5 ng/μl for guts, and 457.9 ± 29.7 ng/μl for dsRNA experimental samples. The OD260/280 ratio of all samples was 1.9–2.1. Single-stranded cDNA was synthesized for each biological sample from 1.0 μg of total RNA using the M-MLV reverse transcription kit (Invitrogen, Carlsbad, CA, USA) and a random N primer (NNNNNN). The cDNA was diluted 10X for the subsequent RT-qPCR studies.

 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. Pair-wise 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 arthropod-active 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 H2O and quantified with a NanoDrop 2000c spectrophotometer and then stored at -20°C.

 Gene Cloning and Primer Design

This study assessed eight reference genes previously commonly used in RT-qPCR analyses and that have been verified as stable reference genes in other species (Table 1). Primers for 18S (AY748147), 28S (DQ202668), 16S (JX896437), and NADH (JQ321839) were designed according to the sequences downloaded from NCBI. Degenerate primers for the remaining four reference genes including Actin, ArgK, EF1A, Tubulin, and one target gene V-ATPase were designed with CODEHOP1 based on conserved amino acid residues among the Coleopteran insect species (Supplementary Table S1). PCR amplifications were performed in 50 μl reactions containing 10 μl 5 × PCR Buffer (Mg2+ Plus), 1 μl dNTP mix (10 mM of each nucleotide), 5 μ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. Amplicons of the expected sizes were purified, cloned into the pCR4-TOPO vector (Invitrogen, Carlsbad, CA, USA), and sent out for sequencing. The confirmation of the reference gene was done by sequence analysis (Supplementary Data Sheet 1), the primers (Table 1) used for RT-qPCR were designed online2 using previously described parameters (Yang et al., 2015).

Primers used for reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR).

Gene Primer sequences (5′–3′) Length (bp) Efficiency (%) R2 Linear regression
18S F:CCAGTAAGCGCGAGTCATAA 103 99.6 0.9975 y = -3.3323x + 13.778
R: GGTCCGAAGACCTCACTAAATC
28S F:TCGAAACGACCTCAACCTATTC 93 113.0 0.9869 y = -3.0462x + 12.955
R: TTGGCACTCTGACCGAAATC
16S F: GGACCTGCCCACTGAATTATTA 106 97.6 0.9992 y = -3.3805x + 15.749
R: TTCTCATCAAACCATTCATACAAGC
EF1A F: CCTGAAGTGGAAGACGAAGAG 98 106.1 0.9782 y = -3.1839x + 20.008
R: AGAAGGAAAGGCTGATGGTAAA
NADH F: AGTAAGAGGAGTAAAGGCATGAAA 75 91.2 0.9961 y = -3.5523x + 20.036
R:CCTATATGGTTGATTGATGATAAGGC
Actin F: GCGTAACCTTCGTAGATTGGTA 100 102.5 0.9994 y = -3.2646x + 19.117
R: CAAGCTGTACTCTCCCTGTATG
Tubulin F: ACAGGTTTCAAAGTGGGTATCA 102 101.6 0.9994 y = -3.2836x + 23.82
R: GGTGGTGTTTGACAACATGC
ArgK F: GTCGAGCTTAGCCTTGTTAGAG 103 110.6 0.9661 y = -3.091x + 23.19
R: GCTGGGTTTCCTCACTTTCT
V-ATPase F: CCTCAAGGTACACCTCCAATTC 84 91.9 0.9935 y = -3.5346x + 20.23
R: AGCTAATGTTCCAGGACGAATAA
Reverse Transcriptase-Quantitative Polymerase Chain Reaction (RT-qPCR)

Details regarding RT-qPCR amplifications and programs were provided in Yang et al. (2014); briefly, PCR reactions (20 μl) contained 7.0 μl of ddH2O, 10.0 μl of 2× SYBR Green Master Mix (Bio-Rad), 1.0 μl of each specific primer (10 μM), and 1.0 μl of first-strand cDNA template. Reactions occurred in 96-well format Microseal PCR plates (Bio-Rad) in triplicate. Reactions were performed in a MyiQ single Color Real-Time PCR Detection System (Bio-Rad). The standard curve for each candidate was generated from cDNA serial dilutions (1/5, 1/25, 1/125, 1/625, and 1/3125). The corresponding RT-qPCR efficiencies (E) were expressed in percentage according to the equation: E = (10[-1/slope]-1) × 100. The coefficients of determination (R2) for a linear regression model with one independent variable was obtained according to the method in excel.

 Determination of Reference Gene Expression Stability

The reference gene expression was evaluated using geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), BestKeeper (Pfaffl et al., 2004), and the ΔCt method (Silver et al., 2006). Candidate reference genes were analyzed with RefFinder3 (Xie et al., 2012). Optimal reference gene number for target gene normalization was determined by pairwise variation (Vn/Vn+1) (0.15 recommended threshold); V-values were computed by geNorm (Vandesompele et al., 2002).

 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.

 Results 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 melt-curve analysis (Supplementary Figure S4). The PCR efficiency (E), correlation coefficient (R2), 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 Cq-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 Cq-values less than 15 cycles. Actin, EF1A, Tubulin, NADH, and V-ATPase had Cq-values ranging from 17 to 24 cycles. 28S and ArgK were the most- and least-expressed reference genes, respectively (Figure 1).

Expression profiles of the eight candidate reference genes and one target gene in all three experiments in Coccinella septempunctata. The whiskers indicate the standard error of the mean. (A–C) indicated the gene expression values under developmental stage, tissue, and dsRNA experimental conditions, respectively

 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.

Stability of candidate reference gene expression under different experimental conditions.

Conditions CRGs∗ geNorm
 Normfider
 BestKeeper
 ΔCt
Stability Rank Stability Rank Stability Rank Stability Rank
Developmental 18S 1.032 6 0.154 1 0.876 6 1.658 4
Stage 28S 0.901 5 0.409 2 0.620 4 1.560 2
16S 0.668 1 0.710 3 0.338 1 1.459 1
EF1A 0.720 2 1.429 6 0.734 5 1.765 6
NADH 0.819 4 0.806 4 0.581 3 1.592 3
Actin 1.440 7 2.209 7 1.893 7 2.697 7
Tubulin 0.743 3 1.145 5 0.578 2 1.670 5
ArgK 1.997 8 3.815 8 3.041 8 3.947 8
Tissue 18S 0.350 1 0.175 1 0.372 3 1.043 4
28S 0.350 1 0.175 1 0.350 2 1.040 3
16S 0.423 2 0.210 2 0.157 1 0.996 1
EF1A 0.597 5 0.830 3 0.845 6 1.287 6
NADH 0.464 3 0.210 2 0.398 5 1.020 2
Actin 0.922 6 1.731 4 1.249 7 1.908 7
Tubulin 0.498 4 0.210 2 0.389 4 1.055 5
ArgK 1.134 7 1.793 5 1.292 8 1.978 8
dsRNA 18S 0.335 1 0.354 5 0.232 2 0.593 5
28S 0.524 7 0.609 8 0.406 8 0.761 8
16S 0.442 5 0.392 6 0.322 6 0.627 6
EF1A 0.393 3 0.200 1 0.258 4 0.563 2
NADH 0.478 6 0.462 7 0.382 7 0.684 7
Actin 0.374 2 0.289 3 0.200 1 0.558 1
Tubulin 0.335 1 0.287 2 0.259 5 0.576 3
ArgK 0.412 4 0.294 4 0.241 3 0.584 4
∗Candidate reference genes.

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 ΔCt 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 ΔCt-based reference gene stability sequences for each experimental condition.

 Comprehensive Ranking of Expression Stability

For the development study, the integrated reference gene rankings (calculated using RefFinder) from most to least stable were as follows: 16S, 28S, NADH, 18S, Tubulin, EF1A, Actin, and ArgK (Figure 2A). In different tissues, the integrated reference gene rankings were: 28S, 16S, 18S, NADH, Tubulin, EF1A, Actin, and ArgK (Figure 2B). For the RNAi treatment, the integrated reference gene rankings were: Actin, Tubulin, EF1A, 18S, ArgK, 16S, NADH, and 28S (Figure 2C).

Stability of candidate reference genes expression in C. septempunctata under different treatment according to their stability value by RefFinder. A lower Geomean value indicates more stable expression. (A–C) indicated the gene stability values under developmental stage, tissue, and dsRNA experimental conditions, respectively.

 Optimal Number of Candidate Reference Genes According to <italic>geNorm</italic>

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).

Pairwise variation (V) values in three groups using geNorm.

 Relative Gene Expression of <italic>V-ATPase</italic>

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 H2O controls (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).

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).

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).

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).

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).

 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 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 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 Ct-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 Cq-values less than 15 cycles, 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, 2016; Vé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.

 Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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

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

Alignment of v-ATPase A ORFs between Coccinella septempunctata and Diabrotica virgifera virgifera. Identical nucleotides are highlighted in black boxes.

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.

The agrose gel profile of the eight candidate reference genes and one target gene in C. septempunctata. M, EZ LoadTM 100 bp Molecular Ruler; Templates in the PCR reactions were as follows: (1) 28S, (2) 18S, (3) 16S, (4) EF1A, (5) Tubulin, (6) V-ATPase, (7) Actin, (8) NADH, and (9) ArgK.

Melting curves of the eight candidate reference genes and one target gene in C. septempunctata.

Standard curves of the eight candidate reference genes and one target gene in C. septempunctata.

 References Alvarez-Alfageme F. Palinkas Z. Bigler F. Romeis J. (2012). Development of an early-tier laboratory bioassay for assessing the impact of orally-active insecticidal compounds on larvae of Coccinella septempunctata (Coleoptera: Coccinellidae). Environ. Entomol. 41 1687–1693. 10.1603/EN12032 Andersen C. L. Jensen J. L. Orntoft T. F. (2004). Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64 5245–5250. 10.1158/0008-5472.CAN-04-0496 Bachman P. M. Bolognesi R. Moar W. J. Mueller G. M. Paradise M. S. Ramaseshadri P. (2013). Characterization of the spectrum of insecticidal activity of a double-stranded RNA with targeted activity against western corn rootworm (Diabrotica virgifera virgifera LeConte). Transgenic Res. 22 1207–1222. 10.1007/s11248-013-9716-5 Baum J. A. Bogaert T. Clinton W. Heck G. R. Feldmann P. Ilagan O. (2007). Control of coleopteran insect pests through RNA interference. Nat. Biotechnol. 25 1322–1326. 10.1038/nbt1359 Berezikov E. (2011). Evolution of microRNA diversity and regulation in animals. Nat. Rev. Genet. 12 846–860. 10.1038/nrg3079 Bustin S. A. Benes V. Garson J. Hellemans J. Huggett J. Kubista M. (2013). The need for transparency and good practices in the qPCR literature. Nat. Methods 10 1063–1067. 10.1038/nmeth.2697 Bustin S. A. Benes V. Garson J. A. Hellemans J. Huggett J. Kubista M. (2009). The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55 611–622. 10.1373/clinchem.2008.112797 Bustin S. A. Benes V. Nolan T. Pfaffl M. W. (2005). Quantitative real-time RT-PCR–a perspective. J. Mol. Endocrinol. 34 597–601. 10.1677/jme.1.01755 Harwood J. D. Samson R. A. Obrycki J. J. (2007). Temporal detection of Cry1Ab-endotoxins in coccinellid predators from fields of Bacillus thuringiensis corn. Bull. Entomol. Res. 97 643–648. 10.1017/S000748530700524X Harwood J. D. Wallin W. G. Obrycki J. J. (2005). Uptake of Bt endotoxins by nontarget herbivores and higher order arthropod predators: molecular evidence from a transgenic corn agroecosystem. Mol. Ecol. 14 2815–2823. 10.1111/j.1365-294X.2005.02611.x Hellemans J. Vandesompele J. (2014). Selection of reliable reference genes for RT-qPCR analysis. Methods Mol. Biol. 1160 19–26. 10.1007/978-1-4939-0733-5_3 Kalushkov P. Hodek I. (2004). The effects of thirteen species of aphids on some life history parameters of the ladybird Coccinella septempunctata. Biol. Control 49 21–32. Kupferschmidt K. (2013). A lethal dose of RNA. Science 341 732–733. 10.1126/science.341.6147.732 Li R. Xie W. Wang S. Wu Q. Yang N. Yang X. (2013). Reference gene selection for qRT-PCR analysis in the sweetpotato whitefly, Bemisia tabaci (Hemiptera: Aleyrodidae). PLoS ONE 8:e53006. 10.1371/journal.pone.0053006 Livak K. J. Schmittgen T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25 402–408. 10.1006/meth.2001.1262 Lundgren J. G. Duan J. J. (2013). RNAi-based insecticidal crops: potential effects on nontarget species. BioScience 63 657–665. 10.1525/bio.2013.63.8.8 Mao Y. B. Cai W. J. Wang J. W. Hong G. J. Tao X. Y. Wang L. J. (2007). Silencing a cotton bollworm P450 monooxygenase gene by plant-mediated RNAi impairs larval tolerance of gossypol. Nat. Biotechnol. 25 1307–1313. 10.1038/nbt1352 Palli S. R. (2014). RNA interference in Colorado potato beetle: steps toward development of dsRNA as a commercial insecticide. Curr. Opin. Insect Sci. 6 1–8. 10.1016/j.cois.2014.09.011 Pan H. Xu L. Noland J. E. Li H. Siegfried B. D. Zhou X. (2016). Assessment of potential risks of dietary RNAi to a soil micro-arthropod, Sinella curviseta Brook (Collembola: Entomobryidae). Front. Plant Sci. 7:1028. 10.3389/fpls.2016.01028 Pan H. Yang X. Bidne K. Hellmich R. L. Siegfried B. D. Zhou X. (2015a). Selection of reference genes for RT-qPCR analysis in the monarch butterfly, Danaus plexippus (L.), a migrating bio-indicator. PloS ONE 10:e0129482. 10.1371/journal.pone.0129482 Pan H. Yang X. Siegfried B. D. Zhou X. (2015b). A comprehensive selection of reference genes for RT-qPCR analysis in a predatory lady beetle, Hippodamia convergens (Coleoptera: Coccinellidae). PLoS ONE 10:e0125868. 10.1371/journal.pone.0125868 Petrick J. S. Frierdich G. E. Carleton S. M. Kessenich C. R. Silvanovich A. Zhang Y. (2016). Corn rootworm-active RNA DvSnf7: repeat dose oral toxicology assessment in support of human and mammalian safety. Regul. Toxicol. Pharmacol. 81 57–68. 10.1016/j.yrtph.2016.07.009 Pfaffl M. W. Tichopad A. Prgomet C. Neuvians T. P. (2004). Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestKeeper–excel-based tool using pair-wise correlations. Biotechnol. Lett. 26 509–515. 10.1023/B:BILE.0000019559.84305.47 Roberts A. F. Devos Y. Lemgo G. N. Y. Zhou X. (2015). Biosafety research for non-target organism risk assessment of RNAi-based GE plants. Front. Plant Sci. 6:958. 10.3389/fpls.2015.00958 Romeis J. Bartsch D. Bigler F. Candolfi M. P. Gielkens M. M. C. Hartley S. E. (2008). Assessment of risk of insect-resistant transgenic crops to nontarget arthropods. Nat. Biotechnol. 26 203–208. 10.1038/nbt1381 Silver N. Best S. Jiang J. Thein S. L. (2006). Selection of housekeeping genes for gene expression studies in human reticulocytes using real-time PCR. BMC Mol. Biol. 7:33–33. 10.1186/1471-2199-7-33 Sinha D. K. Smith C. M. (2014). Selection of reference genes for expression analysis in Diuraphis noxia (Hemiptera: Aphididae) fed on resistant and susceptible wheat plants. Sci. Rep. 4:5059. 10.1038/srep05059 Strube C. Buschbaum S. Wolken S. Schnieder T. (2008). Evaluation of reference genes for quantitative real-time PCR to investigate protein disulfide isomerase transcription pattern in the bovine lungworm Dictyocaulus viviparus. Gene 425 36–43. 10.1016/j.gene.2008.08.001 Tan J. Levine S. L. Bachman P. M. Jensen P. D. Mueller G. M. Uffman J. P. (2016). No impact of DvSnf7 RNA on honey bee (Apis mellifera L.) adults and larvae in dietary feeding tests. Environ. Toxicol. Chem. 35 287–294. 10.1002/etc.3075 Vandesompele J. De Preter K. Pattyn F. Poppe B. Van Roy N. De Paepe A. (2002). Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3:RESEARCH0034. 10.1186/gb-2002-3-7-research0034 Vélez A. M. Jurzenski J. Matz N. Zhou X. Wang H. Ellis M. (2016). Developing an in vivo toxicity assay for RNAi risk assessment in honey bees, Apis mellifera L. Chemosphere 14 1083–1090. 10.1016/j.chemosphere.2015.09.068 Wang H. Chen S. Jiang J. Zhang F. Chen F. (2015). Reference gene selection for cross-species and cross-ploidy level comparisons in Chrysanthemum spp. Sci. Rep. 5:8094. 10.1038/srep08094 Xie F. Xiao P. Chen D. Xu L. Zhang B. (2012). miRDeepFinder: a miRNA analysis tool for deep sequencing of plant small RNAs. Plant Mol. Biol. 80 75–84. 10.1007/s11103-012-9885-2 Xu L. H. Zeng B. S. Norland J. E. Huang Y. P. Zhou X. G. (2015). The coming of RNA-based pest controls. J. Plant Prot. 42 673–690. Yang C. Pan H. Liu Y. Zhou X. (2014). Selection of reference genes for expression analysis using quantitative real-time PCR in the pea aphid, Acyrthosiphon pisum (Harris) (Hemiptera, Aphididae). PLoS ONE 9:e110454. 10.1371/journal.pone.0110454 Yang C. Pan H. Noland J. E. Zhang D. Zhang Z. Liu Y. (2015). Selection of reference genes for RT-qPCR analysis in a predatory biological control agent, Coleomegilla maculata (Coleoptera: Coccinellidae). Sci. Rep. 5:18201. 10.1038/srep18201 Zha W. Peng X. Chen R. Du B. Zhu L. He G. (2011). Knockdown of midgut genes by dsRNA-transgenic plant-mediated RNA interference in the hemipteran insect Nilaparvata lugens. PloS ONE 6:e20504. 10.1371/journal.pone.0020504 Zhu X. Yuan M. Shakeel M. Zhang Y. Wang S. Wang X. (2014). Selection and evaluation of reference genes for expression analysis using qRT-PCR in the beet armyworm Spodoptera exigua (Hubner) (Lepidoptera: Noctuidae). PloS ONE 9:e84730. 10.1371/journal.pone.0084730

http://blocks.fhcrc.org/codehop.html

https://www.idtdna.com/Primerquest/Home/Index

http://fulxie.0fees.us/?ckattempt=1