Selection of Reference Genes for RT-qPCR Analysis in the Monarch Butterfly, Danaus plexippus (L.), a Migrating Bio-Indicator

Reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) is a powerful technique to quantify gene expression. To facilitate gene expression study and obtain accurate results, normalization relative to stably expressed reference genes is crucial. The monarch butterfly, Danaus plexippus (L.), is one of the most recognized insect species for its spectacular annual migration across North America. Besides its great voyages, D. plexippus has drawn attention to its role as a bio-indicator, ranging from genetically modified organisms (GMOs) to natural ecosystems. In this study, nine reference genes from D. plexippus genome were selected as the candidate reference genes. The expression profiles of these candidates under various biotic and abiotic conditions were evaluated using the four readily available computational programs, BestKeeper, Normfinder, geNorm, and ΔCt method, respectively. Moreover, RefFinder, a web-based computational platform integrating the four above mentioned algorisms, provided a comprehensive ranking of the stability of these reference genes. As a result, a suite of reference genes were recommended for each experimental condition. Specifically, elongation factor 1α (EF1A) and ribosomal protein 49 (RP49) were the most stable reference genes, respectively, under biotic (development, tissue, and sex) and abiotic (photoperiod, temperature, and dietary RNAi) conditions. With the recent release of a 273-million base pair draft genome, results from this study allow us to establish a standardized RT-qPCR analysis and lay a foundation for the subsequent genomic and functional genomic research in D. plexippus, a major bio-indicator and an emerging model for migratory animals.


Introduction
The advent of next-generation sequencing technologies has led to a significant increase in transcriptomic and genomic output for various organisms [1][2][3]. Validation of gene expression has become a standard for reporting and assessing the quality of these transcriptomic and genomic resources [4,5]. Reverse transcriptase-quantitative polymerase chain reaction (RT-qPCR) is a powerful technique to target and quantify gene expression [6], however, there remain limitations can significantly influence the normalization of gene expression, such as variations in RNA extraction, RNA quality and integrity, cDNA quality, and PCR efficiency [7][8][9][10]. In RT-qPCR, a commonly used technique to normalize the gene expression data is to introduce one or multiple reference genes, which are stably expressed across various experimental conditions and serve as the internal control [6,9,10]. Although these reference genes have been defined functionally as 'constitutively expressed to maintain cellular function', they do not necessarily meet prerequisites for a good reference gene that can be 'expressed at constant levels across various biotic and abiotic conditions [4,8,10]. Many studies have proved that some commonly used reference genes express differentially across different experimental conditions; researchers have documented that multiple reference genes should be used for accurate normalization [4,8,9,11,12].
RNA interference (RNAi) is a biological process in which RNA molecules inhibit gene expression, typically by causing the destruction of specific mRNA molecules. In view of its high sequence specific reality, RNAi has drawn great attention from the crop protectors who recognized it as a novel and environmentally friendly ways for insect pest control [13,14]. Recently, RNAi-based transgenic crops targeting insects have been developed [15][16][17]. One of the major ecological concerns with respect to the bio-safety of transgenic crops on the environment is their potential effects on non-target organisms [14]. The monarch butterfly is recognized as an important symbol of biodiversity and had been selected as a surrogate species for the risk assessment of Bacillus thuringiensis (Bt) and RNAi-based transgenic crops. Given the nature of RNAi mechanisms, non-target effects could occur via the modulation of gene expressions in these organisms [13]. Therefore, RT-qPCR will be a major research tool to evaluate potential non-target effects of this new biotechnology.
The monarch butterfly, Danaus plexippus(L.) (Lepidoptera: Danainae), is famous for its late summer/autumn southward migration from the United States and southern Canada to Mexico, and northward return in spring, which occurs over the lifespan of three to four generations of the butterfly [18]. In addition to this migration, D. plexippus has drawn attention as a surrogate species for the ecological risk assessment of genetically modified organisms (GMOs) [19][20][21][22] and as an eco-indicator for global climate change [23,24]. The recent release of a 273-million base pair draft genome, including 16,866 protein-coding genes [25], provides an unprecedented opportunity to investigate the genetic basis governing monarch migration, including but not limited to genes involved in their migratory behavior, circadian rhythm, juvenile hormone modulation, and warning coloration [25][26][27]. Recent studies have identified the specific areas in the genome of the monarch that regulate migration. No genetic difference was observed between migrating and non-migrating monarchs but certain genes are expressed solely in migrating monarchs [25,27].
The objective of this study was to determine suitable reference genes with stable expression in D. plexippus across various biotic and abiotic conditions. Here, nine reference genes from D. plexippus genome were selected as candidate reference genes [25], including elongation factor 1α (EF1A), glyceralde hyde-3-phosphate dehydro-genase (GAPDH), nicotinamide adenine dinucleotide (NADH), cyclophilins A (CypA), ribosomal protein S5 (RPS5), ribosomal protein 49 (RP49), vacuolar-type H + -ATPase (v-ATPase), 28S ribosomal RNA (28S), and 18S ribosomal RNA (18S). All of these candidates have been used previously as endogenous references for gene expression analyses [4,[8][9][10][11][12]. Stability of these genes was investigated under three biotic (developmental stage, tissue type, and sex) and three abiotic (temperature, photoperiod, and dietary RNAi) conditions. As a result, different sets of reference genes were recommended for each experimental condition. To validate the recommendations, effectiveness of these candidates were further examined by RT-qPCR analysis of a circadian clock gene timeless (tim) [28].

Ethics Statement
Eggs of Danaus plexippus were collected from common milkweed, Asclepias syriaca, near Ames, Iowa spring and summer 2013. No specific permit was required for the field collection. Danaus plexippus is a common American butterfly species in the United States. The permit to move live plant pests, noxious weeds, and soil were authorized by the United States Department of Agriculture Animal and Plant Health Inspection Service (Permit number: P526P-13-03521). Four to five generations of monarch adults were screened for presence of the protozoan parasite, Ophryocystis elektroscirrha, in order to eliminate it from the population. D. plexippus larvae from the maintenance colony were reared on either common milkweed or tropical milkweed, Asclepias curassavica. Larvae during the experiments were fed honeyvine milkweed, Cynanchum laeve (syn. Ampelamus albidus) at 25 ± 1°C (16L: 8D), and adults were fed 15% sugar solution.

Experimental conditions
Biotic factors. The different developmental stages including eggs, 1 st , 2 nd , 3 rd , 4 th instar larvae (collected at the first day of each instar), pupae, and adults. Tissues, including head, midgut, and carcass, were dissected from the 4 th instar larvae. Male and female adults were collected separately to examine expression profiles of candidate genes between sexes.
Abiotic factors. To examine temperature influence, 2 nd instar D. plexippus were exposed to 10°C, 22°C, and 37°C for 2 h. For photoperiod, 2 nd instars were treated with exposure to 16:8 h, 12:12 h, and 8:16 h light: dark regime for 3 d. For the double-stranded RNA (dsRNA) feeding bioassay, there were three treatments including dsRNAs synthesized from D. plexippus vATPase subunit A (dsDP), β-glucuronidase (GUS) (dsGUS), and water control. Neonates (< 12 h) were given 1.6 μl of water solutions containing 5.0 μg/μl dsRNA each day for two days. Larvae were collected on the third day of the experiment. Altogether, 16μg of dsRNA was ingested by each D. plexippus larva.
For egg samples, six individuals were collected for each replicate; for pupa samples, one pupa was collected for each replicate. For the other biotic and abiotic conditions, approximately five individuals were collected for each treatment with three replicates. All samples were snap frozen with liquid nitrogen in 1.5 ml centrifuge tubes and then stored at -80°C.

Total RNA extraction and cDNA synthesis
Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA) according to previously described methods [29,30]. First-strand cDNA was synthesized from 700 ng of total RNA using the M-MLV reverse transcription kit (Invitrogen, Carlsbad, CA) according to manufacturer's recommendations. The cDNA was diluted 10-fold for the subsequent RT-qPCR studies.

Reference gene primer design and RT-qPCR
Nine reference genes were selected ( Table 1). The primers for the RT-qPCR were designed online (https://www.idtdna.com/Primerquest/Home/Index). The information of RT-qPCR amplification reaction and program were described in detail in our previous study [29,30]. The standard curve and PCR efficiency of each candidate reference gene were constructed and calculated according to previously described methods [29,30].

Data analysis
Stability of the nine candidate reference genes were evaluated by four computational programs, including geNorm [6], NormFinder [10], BestKeeper [31], and the ΔC t method [32]. In addition, RefFinder (http://www.leonxie.com/referencegene.php), a web-based comprehensive platform, which integrates the four above mentioned algorithms, provided the overall ranking of the stability of these reference genes.

Validation of selected reference genes
The photoperiod-relevant tim gene of D. plexippus was used to evaluate the candidate reference genes. tim expression levels were determined in 2 nd instars treated for 3 d under three photoperiod conditions (16L:8D, 12L:12D, and 8L:12D). Relative quantification of tim in different samples was performed using the 2 -ΔΔCt method [33].

Performance and Ct values of candidate reference genes
Each of the tested genes displayed a single amplicon with the expected size on a 1.5% agarose gel (S1 Fig). Moreover, gene-specific of these candidates was affirmed by a single peak in  Table 1. Gene expression analyses of the nine reference genes exhibited a broad C t range (Fig 1). C t values ranged from 7 to 29, while most of the values were distributed between 11 and 22. The most abundant transcripts were 18S and 28S; and the least abundant transcripts were GAPDH and CypA.

Stability of candidate reference genes under biotic conditions
Developmental stage. All software programs, except BestKeeper, identified EF1A as the most stable gene (Table 2).
According to RefFinder, the overall order from the most stable to the least stable reference genes across different developmental stages was: EF1A, 28S, RPS5, RP49, v-ATPase, 18S, NADH, GAPDH, CypA (Fig 2A). geNorm analysis revealed that the pair-wise variation value V3/4 was 0.154, close to the proposed 0.15 cut-off (Fig 3), suggesting that three reference genes were required for reliable normalization throughout developmental stages.
Tissue. Majority of the software programs ranked 28S as the most stable gene (Table 2). Based on RefFinder, the overall order from the most to the least stable reference genes in different tissue was: 28S, EF1A, CypA, NADH, RP49, v-ATPase, RPS5, GAPDH, 18S (Fig 2B). geNorm analysis revealed that the first V-value < 0.15 showed at V3/4 (Fig 3), suggesting that three reference genes were sufficient for reliable normalization in different tissue types.
Sex. According to the results generated by geNorm, 18S was ranked as the most stable gene whereas the other three programs identified it as the second stable gene ( Table 2). According to RefFinder, the overall order from the most stable to the least stable reference genes in both female and male adults was: 18S, RP49, 28S, v-ATPase, EF1A, GAPDH, RPS5, CypA, NADH (Fig 2C). geNorm analysis revealed that the first V-value < 0.15 showed at V2/3 (Fig 3), suggesting that two reference genes were enough for accurate normalization in female and male adults.

Stability of candidate reference genes under abiotic conditions
Temperature. All software programs identified EF1A as the most stable gene (Table 3). According to RefFinder, the overall order from the most stable to the least stable reference genes under the temperature stress was: EF1A, RP49, 18S, v-ATPase, CypA, RPS5, 28S, NADH, GAPDH (Fig 2D). geNorm analysis revealed that first V-value < 0.15 showed at V2/3 (Fig 3), suggesting that two reference genes were sufficient for accurate normalization across different temperatures.
Photoperiod. All software programs identified EF1A as the most or second to most stable gene (Table 3). According to RefFinder, the overall order from the most stable to the least stable reference genes under the photoperiod stress was: EF1A, v-ATPase, 18S, RP49, CypA, RPS5, 28S, NADH, GAPDH (Fig 2E). geNorm analysis revealed that first V-value < 0.15 showed at Determination of the optimal number of reference genes. Pairwise variation (V) 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 [6].
dietary RNAi. All software programs, except BestKeeper, identified GAPDH as the most stable gene (Table 3). According to RefFinder, the overall order from the most stable to the least stable reference genes under the dietary RNAi treatment was: GAPDH, EF1A, 28S, 18S, RP49, RPS5, NADH, CypA, v-ATPase (Fig 2F). geNorm analysis revealed that the pair-wise variation value V4/5 was below the proposed 0.15 cut-off (Fig 3). This result suggests that four reference genes were required for accurate normalization under the dietary RNAi condition.

Validation of selected reference genes
Using two best reference gene combinations for normalization, two genes [EF1A and v-ATPase; NF (1-2)] or three genes [EF1A, v-ATPase, and 18S; NF (1-3)], similar expression levels of tim occurred within each photoperiod (Fig 4). Under the 12L: 12D photoperiod condition, the expression level of tim was the highest when the reference gene with the highest Geomean value (GAPDH; NF9) was used (Fig 4).

Discussion
The monarch butterfly is commonly used as a non-target insect for the ecological risk assessment of transgenic crops in the U.S.; plus interest in the migratory biology and genetics of this insect suggest future post-genome studies will occur. In this study, expression profiles of nine candidate reference genes from D. plexippus were evaluated under diverse experimental conditions. This research is the first step toward establishing a standardized RT-qPCR analysis for this insect.
Our results demonstrate that the best-suited reference genes can vary depending on biotic and abiotic factors (Tables 2 and 3). Thus, reference genes need to be selected based on the , and NF9 indicate that the expression of tim was normalized using the best, the top two, the top three, or, the worst reference genes, respectively. Bar represents the mean and standard error of three biological replicates.
doi:10.1371/journal.pone.0129482.g004 experimental conditions. Unfortunately, a universal reference gene that covers all experimental conditions is unrealistic. For example, EF1A, which plays an important role in translation by catalyzing the GTP-dependent binding of aminoacyl-tRNA to the acceptor site of the ribosome, is the most stable reference gene for biotic factors (Fig 2G). This result is consistent with previous studies in other insects [34,35]. However, EF1A ranked poorly (7 th ) as a reference gene for abiotic factors (Fig 2H).
To avoid biased normalization, many researchers have started to advocate the use of multiple reference genes to analyze gene expression [4,5,9,11,12,29,30,36]. Our results demonstrate that two reference genes are required for reliable normalization under photoperiod, temperature, and sex conditions; three reference genes are required for reliable normalization in different tissue types and under different developmental stages; four reference genes are sufficient to normalize the expression and provide a more conservative estimation of target gene expression under dietary RNAi condition (Fig 3). As a result, we suggest that using different combinations of reference genes are necessary for studying gene expression in D. plexippus, especially with variable experimental conditions. This work represents an important first step toward establishing a standardized gene analysis framework for D. plexippus. Analysis of our data indicate that EF1A and RP49 are the most stable reference genes, respectively, under biotic (development, tissue, and sex) and abiotic (photoperiod, temperature, and dietary RNAi) conditions. With the recent release of D. plexippus genome, results from this study provide a critical piece for the subsequent genomics and functional genomics research in this emerging insect model, and shed light on the ecological risk assessment of RNAi-based biopesticides on this non-target bio-indicator agent.