Reference gene selection for RT-qPCR analysis in Harmonia axyridis, a global invasive lady beetle

Harmonia axyridis is a voracious predator, a biological control agent, and one of the world most invasive insect species. The advent of next-generation sequencing platforms has propelled entomological research into the genomics and post-genomics era. Real-time quantitative PCR (RT-qPCR), a primary tool for gene expression analysis, is a core technique governs the genomic research. The selection of internal reference genes, however, can significantly impact the interpretation of RT-qPCR results. The overall goal of this study is to identify the reference genes in the highly invasive H. axyridis. Our central hypothesis is that the suitable reference genes for RT-qPCR analysis can be selected from housekeeping genes. To test this hypothesis, the stability of nine housekeeping genes, including 18S, 28S, ACTB, ATP1A1, GAPDH, HSP70, HSP90, RP49, and ATP6V1A, were investigated under both biotic (developmental time, tissue and sex), and abiotic (temperature, photoperiod, in vivo RNAi) conditions. Gene expression profiles were analyzed by geNorm, Normfinder, BestKeeper, and the ΔCt method. Our combined results recommend a specific set of reference genes for each experimental condition. With the recent influx of genomic information for H. axyridis, this study lays the foundation for an in-depth omics dissection of biological invasion in this emerging model.

The multicolored Asian lady beetle, Harmonia axyridis (Coleoptera: Coccinellidae), a generalist predator, preys on aphids and scale insects on crops and other plants 1 . Harmonia axyridis is native to central and eastern Asian. To exploit its ecosystem services, numerous releases were attempted in North America and Europe, as early as 1916 2,3 . Due to its broad range of preys and incredible consumption rate, H. axyridis indeed has been used to control aphids [4][5][6] and other sap-sucking arthropod pests 7,8 . However, the worldwide propagation of H. axyridis threatens the indigenous lady beetles and other non-target species [9][10][11] . Considered as "the most invasive ladybird on Earth", the role of H. axyridis has shifted from a global biological control agent to an invasive alien species 12 . Multiple factors contribute to this transition. Predation of eggs and larvae of other lady beetle species is one of the reasons which leads to the decline of native species 13,14 . A higher level of resistance to infection is the other major reason to benefit its competition in the field [15][16][17] . The molecular basis of this resistance, however, is poorly understood.
Double-stranded RNA (dsRNA) can induce sequence-specific posttranscriptional gene silencing in many organisms, i.e., RNA interference (RNAi) 18,19 . RNAi can not only investigate gene functions in vivo or in vitro, but also offers a novel approach with a brand new mode-of-action to control arthropod pests [20][21][22][23][24] . With a recent influx of genomic information for H. axyridis, there is an increasing need for the development of genetic tools to functionally interpret the sequencing data 20,[24][25][26] .
Real-time quantitative PCR (RT-qPCR) has been used primarily for gene expression quantification [27][28][29] . RT-qPCR analysis is highly sensitive, and its accuracy can be affected by RNA quantity, transcription efficiency, amplification efficiency and experimental procedures between samples. To avoid biases, normalization of gene expression is an essential step 30 . The most common practice is to compare a target gene expression with an internal reference gene 31 . Housekeeping genes, such as beta-actin (ACTB), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), and translation elongation factor 1-alpha (EF1A) 32,33 have been used extensively for RT-qPCR analysis.
However, under any given experimental condition, the expression of these commonly used reference genes may vary substantially [34][35][36][37] . A systematic and customized study for each tested species is recommended for identifying appropriate reference genes 38,39 .
The overall goal of this study is to identify candidate reference genes in the highly invasive H. axyridis. Our objective is to determine the suitable reference genes for RT-qPCR analysis in H. axyridis from selected housekeeping genes, an array of constitutively expressed genes maintaining the basic cellular functions in an organism. We evaluated the stability of nine housekeeping genes under selected biotic and abiotic conditions, respectively. The candidate genes include 18S ribosomal RNA(18S), 28S ribosomal RNA (28S), Na+/K+-ATPase subunit alpha 1 (ATP1A1), heat shock protein 70 (HSP70), heat shock protein 90 (HSP90), ribosomal protein 49 (RP49), V-ATPase subunit A (ATP6V1A), ACTB and GAPDH from H. axyridis. All these housekeeping genes have been used empirically as the reference genes for RT-qPCR analyses in other organisms, especially in insects. The specific environmental conditions range from biotic (developmental stage, tissue type, and sex) to abiotic treatments (temperature, photoperiod, and in vivo RNAi). As a result, a specific set of reference genes is recommended for each given condition.

Results
RT-qPCR analysis. For each candidate reference gene, a single amplicon was produced, as detected by agarose gel electrophoresis analysis and the melting curve analysis. Nonspecific bands were not found, and a single peak was observed in the melting curve analysis. A standard curve was generated for each gene, using a five-fold serial dilution of the pooled cDNA. Efficiency of RT-qPCR ranged between 90 and 110% (Table 1), which is considered standard 40 . Ct values of the nine candidate reference genes ranged from 8 to 27, covering all the experimental conditions (Fig. 1). While the vast majority of Ct values were found between 17 and 26, 18S was the most abundant transcript. ATP1A1, VATP6V1A, and RP49 were the least abundant candidate reference genes.
Recommended reference genes. For repeatable and consistent results, multiple normalizers (≥2 reference genes) are required for RT-qPCR analysis. GeNorm analysis evaluated all pairwise variations under each experimental conditions (Fig. 3). According to Vandesompele et al. 31 , a Vn/Vn + 1 cutoff value of 0.15 means the addition of n + 1 reference gene is not necessary, i.e., the first n references genes are sufficient to normalize qRT-PCR results. The optimal number of reference genes was recommended in Tables 2 and 3, respectively, for biotic and abiotic conditions. Specifically, for different developmental stages, the recommended reference genes were 18S, HSP70, and 28S. For different tissues, the recommendation was 28S, 18S, and RP49. For different sexes, the recommendation was HSP90 and RP49. For different temperature treatments, the recommendation was 18S, 28S, and GAPDH. For different photoperiods, the recommendation was 18S, 28S, and HSP90. Finally, for in vivo RNAi, the best combination was RP49 and ATP1A1.
Validation of selected reference genes. The expression of TPS, a target gene, was evaluated to validate the recommended reference genes under different temperature treatments. Using the most stable reference gene 18S (NF 1), the top two stable reference genes 18S and 28S (NF 1-2), and the top three stable reference genes, 18S, 28S, and GAPDH (NF 1-3) for normalization, TPS expression profiles were similar throughout all three temperature regimes (Fig. 4). In comparison, when ATP6V1A, the least stable candidate (NF 9), was used as the reference gene, TPS expression patterns were inconsistent across different temperature treatments. Specifically, TPS expression was numerically higher at 10 °C, and lower at 22 and 30 °C (Fig. 4).

Discussion
RT-qPCR has been used extensively for quantification of mRNA expression and is a primary tool for genetic research. Although multiple factors, such as RNA extraction, storage, cDNA synthesis, and handling of materials and reagents, can affect the RT-qPCR analysis, a reliable reference gene (set) to overcome confounding variations in an empirical dataset is of particular importance. Normalization by internal controls is an integral part of the quantification process. A single or multiple stably expressed reference genes are required for the normalization process to achieve accurate and reliable results. Each candidate reference gene should be evaluated under specific experimental conditions to ensure a constant level of expression 35 . Following the "Minimum Information for Publication of Quantitative Real-Time PCR Experiments" (MIQE) guideline 41 , reference gene selection study has been carried out for many insect species 34,42,43 , and has become a routine practice to standardize RT-qPCR analysis.
Due to different algorithms, stability ranking derived from the four analytical tools can vary. For example, when H. axyridis was injected with dsRNAs (in vivo RNAi), 28S was rated as the best reference gene by BestKeeper, RP49 was considered as the most stable by Normfinder as well as ΔCT method, whereas ATP1A1 and GAPDH were the top choice by geNorm. Despite some discrepancies in individual rankings, RP49 and ATP1A1 were consistently exhibited a higher level of stability than the rest of the candidates projected by all four algorithms (Table 3), suggesting the importance of (1) using a comprehensive analysis to interpret the dataset and (2) adopting the multiple instead of a single normalizer for RT-qPCR analysis.
In recent years, researchers have been more receptive to use multiple reference genes to replace a single normalizer in RT-qPCR analysis 44 . The optimal number of reference genes is typically determined by geNorm. In this study, three reference genes for recommended for different developmental stages (18S, HSP70, and 28S), tissues (28S, 18S, and RP49), temperatures (18S, 28S and GAPDH), and photoperiods (18S, 28S and HSP90), while two reference genes were required for the reliable normalization in different sexes (HSP90 and RP49), and in vivo RNAi (RP49 and ATP1A1). Our combined results are, in part, consistent with previous studies of other Coccinellidae predatory species (Table 4), especially for ribosome RNAs (rRNAs).
Not surprisingly, rRNAs (e.g., 18S and 28S) were consistently stably expressed throughout the vast majority of biotic and abiotic conditions among the four Coccinellidae species, including H. axyridis, Hippodamia convergens 45 , Coleomegilla maculate 46 , and Coccinella septempunctata 47 . The over-representation of rRNAs in the total RNA pool (>80%), however, can potentially mask the subtle changes of the target gene expression 48 . Therefore, customized reference gene study is still a prerequisite for standardized RT-qPCR analysis in predatory lady beetles. A large body of works has demonstrated that there are no "universal" reference genes applicable for all cell and tissue types and various experimental conditions 49 . As a major structural protein, Actin has been used extensively as the internal control without any validation. In this study, however, Actin was one of the least stable candidates under both biotic and abiotic conditions, except the temperature treatment, which is consistent with the other three Coccinellidae species [45][46][47] .
This study not only provides a standardized procedure for the quantification of gene expression, but also lays a foundation for the genomics and functional genomics dissection of H. axyridis, an emerging model in invasion biology 50 .

Materials and Methods
Insects. Harmonia axyridis was originally collected from the University of Kentucky North Farm (38°07′N, 84°30′W). Harmonia axyridis colony was maintained at 23 ± 1 °C, 12 L:12D photoperiod, 50% relative humidity, and provisioned with pea aphids and sugar water for more than two months. Pea aphid clones were a gift from Dr. John Obrycki (University of Kentucky) and maintained on seedlings of fava beans in a glasshouse. Experimental conditions. Biotic conditions. The developmental stages include eggs (N = 15), four larval instars (N = 5 for each instar, respectively), pupae (N = 1), and adults (one male and one female). Sex of adult beetles was determined by the presence or absence of the male genitalia. Tissues, including head, midgut, and carcass, were dissected from the fourth instar larvae (N = 5).

Biotic Conditions
Abiotic conditions. To examine the effects of temperature, third instars were exposed to 10, 22, and 30 °C for 3 hours. For photoperiod, third-instar larvae were treated with a series of light and dark regime of 16 L:8D, 12 L:12D, and 8 L:16D for two days. For in vivo RNAi, H. axyridis ATP6V1A was the intended molecular target.       . The existence of one peak in melting curve analysis was used to confirm gene-specific amplification and to rule out non-specific amplification and primer-dimer generation. The RT-qPCR was determined for each gene using slope analysis with a linear regression model. Relative standard curves for the transcripts were generated with a serial dilution of cDNA. The corresponding RT-qPCR efficiencies (E) was calculated according to the equation: Validation of selected reference genes. Trehalose-6-phosphate synthase (TPS), the intermediate of trehalose, is a key component in insect energy metabolism and resilience 25,51,55 . The stability of candidate reference genes was evaluated using TPS as the target gene. TPS expression levels under different temperature treatments were calculated based on selected sets of candidate reference genes. Two separate normalization factors (NFs) have been computed based on (1) the geometric mean of the genes with the lowest Geomean values (as determined by RefFinder), and (2) a single normalizer with the lowest or highest Geomean value. Relative expression of TPS in different samples was calculated using the 2 −ΔΔCt method 56 .