Screening for stable internal reference genes for quantitative PCR analysis of Wolbachia-host interactions in whitefly Bemisia tabaci (Homoptera: Aleyrodidae)

Stable reference genes (RGs) determine the reliability of quantitative polymerase chain reaction (qPCR) analyses and it is recommended that different reference genes are used for different types of DNA and tissues. The present study aimed to screen for stable RGs for the qPCR analysis of the immune responses of the whitefl y Bemisia tabaci to the Wolbachia wMel strain from Drosophila melanogaster. A total of eight candidate RGs were evaluated using fi ve different methods, i.e., Coeffi cient of Variation analysis, GeNorm, NormFinder, BestKeeper and ΔCt. The stability of these RGs was assessed for both genomic DNA (gDNA) and complementary DNA (cDNA). The results indicate that β-actin (Actin) and elongation factor 1 alpha (EF-1α) were the most stable RGs for gDNA, whereas 18S rRNA (18S) and glyceraldehyde phosphate dehydrogenase (GAPDH) were the least stable; in contrast, Actin and GAPDH were the most stable for cDNA, whereas RPL29 and ATPase were the least stable. The effectiveness of the most stable RGs was then validated against the least stable using qPCR analysis of the titre of wMel (gDNA) and the transcriptional responses of the antimicrobial peptide Alo-3-like and the phosphatidylinositol-bisphosphate 3-kinase catalytic subunit delta isoform (cDNA) to wMel transfection. The results support the notion that reliable RGs are essential for a qPCR analysis of samples of both gDNA and cDNA. * Corresponding author: e-mail: zxli@cau.edu.cn INTRODUCTION Real-time quantitative PCR (qPCR) is a commonly used technique for gene quantitation at both genomic DNA (gDNA) and complementary DNA (cDNA) levels, which recently has been used in many studies (Artico et al., 2010; Derveaux et al., 2010; Hindson et al., 2013). Nevertheless, the reliability of qPCR is affected by many different factors (Bustin, 2000; Derveaux et al., 2010) and the use of reference genes (RGs) can greatly improve the accuracy of qPCR results by standardizing or normalizing experimental data from different developmental stages, organs or tissues (Artico et al., 2010; Arya et al., 2017). Screening for stable RGs is reported for different species of insects, such as Aphis glycines (Bansal et al., 2012), Nilaparvata lugens (Yuan et al., 2014), Bactericera cockerelli (Ibanez & Tamborindeguy, 2016) and Lipaphis erysimi (Koramutla et al., 2016). The RGs have to be carefully evaluated for particular situations, as it is unlikely that the same RGs (non-evaluated) are generally suitable (Thellin et al., 1999; Vandesompele et al., 2002; Ponton et al., 2011; Zhou & Li, 2016). Several algorithms have been developed for Eur. J. Entomol. 116: 402–412, 2019 doi: 10.14411/eje.2019.041


INTRODUCTION
Real-time quantitative PCR (qPCR) is a commonly used technique for gene quantitation at both genomic DNA (gDNA) and complementary DNA (cDNA) levels, which recently has been used in many studies (Artico et al., 2010;Derveaux et al., 2010;Hindson et al., 2013). Nevertheless, the reliability of qPCR is affected by many different factors (Bustin, 2000;Derveaux et al., 2010) and the use of reference genes (RGs) can greatly improve the accuracy of qPCR results by standardizing or normalizing experimental data from different developmental stages, organs or tissues (Artico et al., 2010;Arya et al., 2017). Screening for stable RGs is reported for different species of insects, such as Aphis glycines (Bansal et al., 2012), Nilaparvata lugens (Yuan et al., 2014), Bactericera cockerelli (Ibanez & Tamborindeguy, 2016) and Lipaphis erysimi (Koramutla et al., 2016). The RGs have to be carefully evaluated for particular situations, as it is unlikely that the same RGs (non-evaluated) are generally suitable (Thellin et al., 1999;Vandesompele et al., 2002;Ponton et al., 2011;Zhou & Li, 2016). Several algorithms have been developed for adenosine triphosphate enzyme (ATPase) ( Table 1). Most of these genes have been used as internal references for real-time qPCR analysis in previous studies (Zhou et al., 1998;Li et al., 2013;Zhou & Li, 2016). These RGs were evaluated against both gDNA and cDNA in the present study, but HSP20 and HSP70 were not used for cDNA samples as it is likely they are induced by Wolbachia transfection.

Quantitative PCR analysis
The specifi city of the primers for the eight RGs (Table 1) was examined using PCR in a total reaction volume of 25 μL containing 2 μL DNA template, 1 μL of each primer (10 μM), 2.5 μL 10 × PCR buffer (TransGen Biotech, Beijing, China), 2 μL dNTPs (2.5 mM), 0.5 μL Taq DNA polymerase (5.0 U/μL) and 16 μL ddH 2 O. The thermocycling program was: 94°C for 5 min, 35 cycles of 94°C for 30 s, 60°C for 30 s, 72°C for 30 s and a fi nal 10 min extension at 72°C. PCR products were analysed using 1% agarose gel electrophoresis. For the qPCR analysis, the DNA templates (gDNA or cDNA) were successively diluted by 5 0 , 5 1 , 5 2 , 5 3 and 5 4 times for construction of standard curves and melting curves. A suitable pair of primers for qPCR analysis should have a correlation coeffi cient (R 2 ) > 0.99, an amplifi cation effi ciency (E) > 90% and a unimodal melting curve (Livak & Schmittgen, 2001). The reaction was performed in a total volume of 20 μL containing 10 μL AceQ ® SYBR ® Green Master Mix (Vazyme, Nanjing, China), 0.4 μL of each primer (10 μM), 1 μL DNA template and 8.2 μL ddH 2 O on ABI 7500 platform (Applied Biosystems, Foster City, California, USA); the primer pair wspQ384/wspQ513 (Table 1) was used to detect the titre of wMel strain in B. tabaci after transfection. The thermocycling program was 50°C for 2 min, 95°C for 5 min, 40 cycles of 95°C for 10 s and 60°C for 30 s. The program used for determining the melting curve was 95°C for 15 s, 60°C for 60 s and 95°C for 15 s. DNase/ RNase-free water was used as the negative control. Each treatment was performed in triplicate.

Stability of RGs
The cycle threshold (Ct) values of RGs for different generations were compared in order to evaluate their stability. The working concentration of gDNA for the qPCR analysis was 20 ng/μL and that of cDNA was 100 ng/μL. The stability of RGs was evaluated using four different algorithms: GeNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), BestKeeper (Pfaffl et al., 2004 and ΔCt method (Silver et al., 2006). The four algorithms rated the stability according to different variables: GeNorm provides an M-value based on the average paired expression ratio. The lower the M-value, the more stable the expression: M ≤ 1.5 was considered to be stable. NormFinder calculates stability based on the inter-and intra-group variance of each candidate RG; the lower the NormFinder value, the more stable the RG. BestKeeper determines the stability of RG according to the standard deviation (SD): The higher the SD, the more unstable the expression. ΔCt method performs pairwise multiple comparisons between the expression levels and identifi es the most stable RGs. Each algorithm sorted RGs according to their stability and assigned the RGs a series of consecutive integers (starting at one) and the geometric average (geomean) of the weights based on the four algorithms was then calculated for each RG. The RG with the lowest geo mean value was considered to be the most stable and the RGs were ranked accordingly. Subsequently, the Coefficient of Variation (CV) was used to measure the population variance of each gene as it is the only method that is not affected by other factors (Boda et al., 2009). Here a threshold of CV = 50% was defi ned: genes with a CV value above this threshold were enous Wolbachia strains were successfully transferred into B. tabaci by microinjection and transinfected isofemale lines were established in the laboratory (Zhong & Li, 2014;Hu & Li, 2015;Zhou & Li, 2016). Recently, we sequenced the transcriptomes of B. tabaci in response to Wolbachia transfection using the RNA-seq technique, which is necessary for a detailed investigation of the host-Wolbachia interaction, including qPCR analysis of the titre of Wolbachia and expression levels of differentially expressed genes (DEGs) and functional analysis of candidate genes via RNA interference (RNAi). All of these studies require the selection of suitable RGs. In a previous study, 15 candidate RGs were evaluated for B. tabaci using the algorithms geNorm and Normfi nder , which revealed that the stability of RGs was greatly affected by both biotic and abiotic factors, and that different RGs should be used depending on the species and conditions. The purpose of the present study was to screen for suitable RGs for qPCR analysis of both gDNA and cDNA from B. tabaci artificially transfected with exogenous Wolbachia. Based on our studies, different RGs are needed for determining the titre of Wolbachia or level of expression of functional genes in different generations when B. tabaci is transfected with the wMel Wolbachia strain from the fruit fl y Drosophila melanogaster.

Insect rearing and transfection through microinjection
The whitefl y B . tabaci (Q cryptic species) was collected in tomato greenhouses of Jinan Academy of Agricultural Sciences, Shandong, China, a nd then maintained on pot-grown plants of cotton in a laboratory (14L : 10D, RH 60-80%, 27 ± 1°C). The 4 th -instar nymphs (pseudopupae) of B. tabaci were artifi cially transinfected with the wMel Wolbachia strain by microinjection as described by Zhou & Li (2016). The wMel strain was isolated from D. melanogaster using the Percoll density-gradient centrifugation method. The primer pair 81F/691R (Table 1) was used to verify the existence of the wMel strain, and on each occasion a volume of 46 nL Wolbachia suspension was injected. After transfection, isofemale lines were established and samples were collected from subsequent generations (G).

Extraction of genomic DNA (gDNA) and synthesis of complementary DNA (cDNA)
Twenty fresh whitefl ies (adults) were used for extraction of gDNA from samples of G 0 , G 2 , G 4 , G 6 and G 7 individuals using the KAC method as described by Zhong & Li (2013) and the purity and concentration of gDNA were checked using a Nan-oDrop ND-2000 (Thermo Scientifi c, Wilmington, DE). Samples of thirty fresh whitefl ies were used for the extraction of total RNA in a 1.5-mL centrifuge tube containing TRIZol reagent (TransGen Biotech, Beijing, China), which were from G 0 , G 1 and G 4 of the transinfected isofemale lines, with the wild type as the control. The fi rst-strand cDNA was synthesized from 500 ng of total RNA using the PrimeScript TM RT reagent Kit (TaKaRa, Beijing, China) according to the supplier's instructions.

Validation of RGs
The reliability of candidate RGs was validated by a qPCR analysis of the titre of wMel (gDNA) and the transcriptional response of two genes (antimicrobial peptide Alo-3-like; pho sphatidylinositol-bisphosphate 3-kinase catalytic subunit delta isoform, Pl3K) to wMel transfection (cDNA) (the primer sequences are included in Table 1). Alo-3-like is an antimicrobial peptide involved in immune response of the host to wMel transfection, and Pl3K is involved in the immune-related signalling pathway in B. tabaci based on our transcriptome sequencing data (Tables S1 and S2; Figs S1 and S2). Here we used the two most stable and two least stable RGs based on their geomean values to normalize the qPCR analysis before the results were compared. The procedures for qPCR analysis were the same as described above.

Data analysis
The statistical differences were analysed using One-way ANOVA followed by Student Newman Keuls (SNK) test at 0.05, 0.01 and 0.001 levels of probability on SPSS v.20.0 (SPSS Inc., Chicago, IL, USA).

Specifi city, amplifi cation effi ciency and cycle thresholds of RG primers
The specifi city of the RG primers (Table 1) was determined using PCR and agarose gel electrophoresis, which indicated that all primer pairs resulted in single specifi c bands of expected size. In addition, the melting curves of all the primer pairs had single peaks, and the correlation coeffi cients and amplifi cation effi ciencies of the standard curves were all within the range of R 2 > 0.99 and 90% < E < 110% (Table 1; Fig. S3). Moreover, the determination of RG using qPCR revealed that the Ct-value of 18S was the lowest (the log of starting template copies) for both gDNA and cDNA templates, while that of GAPDH was the highest for cDNA ( Fig. 1) and that of RPL29 was the most variable (Fig. 1A).

Stability of candidate RGs
Based on qPCR analysis of gDNA and cDNA samples taken from transinfected B. tabaci adults and the different algorithms (GeNorm, NormFinder, BestKeeper or ΔCt method), stabilities were assigned to each RG and the geomeans of these stabilities calculated, which were used to determine their combined ranking ( Table 2). The results indicated that ATPase and EF-1α were the two most stable RGs for gDNA templates and RPL29 and GAPDH the two least stable RGs; in contrast, 18S and Actin were the two most stable RGs for cDNA templates, whereas EF-1α and ATPase were the two least stable RGs.
The CV analysis was the most appropriate method for identifying the most variable RGs, which were excluded from further assessments. The analysis showed that RPL29 and EF-1α were the most variable RGs for both gDNA and cDNA samples (Table 3). After removing RPL29 and EF-1α from the ranking analysis, a corrected stability ranking value was generated for each RG using the different algorithms and fi nally all RGs were ranked according to their geomean values (Table 4). The results indicate that Actin and EF-1α were the two most stable RGs for gDNA samples, whereas Actin and GAPDH were the two most stable for cDNA samples. Interestingly, 18S/GAPDH and 18S/RPL29/ATPase were identifi ed as the least stable RGs for both gDNA and cDNA samples.

Validation of stable RGs
The qPCR analysis of the relative titre (gDNA) of the Wolbachia wMel strain was normalized by using a combination of the two most stable (Actin and EF-1α) and two least stable (18S and GAPDH) RGs. The results revealed that the measurement of Wolbachia titre was indeed affected by the nature of the RG in all of the fi ve generations assessed in this study (Fig. 2). Specifi cally, the relative titres normalized by unstable RGs varied dramatically between generations, especially in G 6 for which the titre was several-fold higher than in that recorded in other generations. In contrast, the titres normalized by stable RGs looked more normal.
The qPCR analysis of the relative expression (cDNA) of two immune-related genes (AMP and Pl3K) was also normalized using a combination of the two most stable (Actin and GAPDH) and two least stable (RPL29/ATPase) RGs. The results indicate that the relative level of expression of AMP (Alo-3-like gene) when normalized by unstable RGs was very high in the wild type, while there were no signifi cant differences between G 0 and G 4 and the control (G 1 ) (Fig. 3A). In the case of Pl3K, the relative level of expression normalized by unstable RGs were exceptionally high in the wild type and in G 0 of the transinfected lines (Fig.  3B). Moreover, the relative levels of expression of both the Alo-3-like gene and Pl3K when normalized using stable RGs were generally consistent with the transcriptome data, whereas those normalized using unstable RGs were not (Tables S1 and S2).

DISCUSSION
Our studies confi rmed that the reference genes signifi cantly affected the accuracy of qPCR analysis of both gDNA and cDNA samples. Moreover, the selection of candidate RGs based on the algorithms used in this study seemed very effective: the scoring and ranking procedures ensured the selection of suitable RGs for the qPCR analysis of both gDNA and cDNA templates. Our results support  the notion that the use of suitable RGs is essential for a qPCR analysis.
The wMel Wolbachia strain from D. melanogaster is able to establish and induce strong cytoplasmic incompatibility (CI) in B. tabaci (Zhou & Li, 2016), an important worldwide agricultural pest, which indicates it might be possible to use it to control this pest. Nevertheless, the titre of wMel fl uctuated during transgenerational transmission in B. tabaci. For instance, the titre of wMel was extremely low after transfection in G 1-2 and G 2 and higher in G 4-6 . This observation prompted us to speculate that the titre of Wolbachia might be modulated by the interaction between the host and Wolbachia. Obviously, clarifying the quantitative relationships between the titre of Wolbachia and the levels of expression of the candidate genes involved in the host-Wolbachia interaction would help us understand the mechanisms underlying the fl uctuations in the Wolbachia titre, especially when the titre is potentially related to the pest control capability of Wolbachia (Breeuwer & Werren, 1993;Noda et al., 2001). Therefore, it is essential to quantitatively measure both the titre of Wolbachia (gDNA) and the levels of expression (cDNA) of functional genes, and thus the screening for stable RGs for both gDNA and cDNA samples is a prerequisite for a reliable qPCR analysis. In the present study, a combination of algorithms was used to evaluate the candidate RGs, which circumvented the drawbacks of using a single alg orithm, and the results indicate that use of the stable RGs rather than unstable RGs could result in more reliable qPCR data. Indeed, the stable RGs specifi cally identifi ed in this study were successfully used in the qPCR analysis of both gDNA and cDNA samples from B. tabaci artifi cially transfected with Wolbachia, which will facilitate the study of the host insect-Wolbachia interaction.
Different algorithms were developed for selecting RGs in the past, but as mentioned above, each algorithm has advantages and disadvantages, which potentially generate biased results. The use of a combination of algorithms may hopefully counteract any bias. For instance, GeNorm and BestKeeper are based on paired comparisons, and the selection of the most appropriate RGs are based on the change in the expression of genes, which does not eliminate the infl uence of co-regulation; on the other hand, Nor-mFinder and the ΔCt method can neutralize this infl uence. Here we used geomeans to aggregate the RGs and generated an integrated ranking for the candidate RGs and the most variable RGs based on the CV analysis were not used in this study. Our corrected results indicate that the combined analysis worked well in selecting stable RGs.
The overall design of our experiments was based on the results of our previous study of the wMel Wolbachia strain after transfection by microinjection, which revealed the titre of wMel differed greatly in the different generations. As noted above, the change in titre was most marked from G 0 to G 6 ; hence, the samples (both gDNA and cDNA) were collected during this period of time. Hopefully, the change in the expression of the functional genes selected and the titre of Wolbachia were effectively detected. In addition, the selection of immune-related genes for validating the RGs was based on those genes that were identifi ed by the transcriptome sequencing, which revealed that an active immune response was induced after infection with Wolbachia, and thus they will be functionally analysed in the next step. The antimicrobial peptide Alo-3 was greatly and EF-1α are excluded from gDNA and cDNA samples, respectively, in the corrected ranking analysis based on their highest CV values (Table 3).
down-regulated, while the PL3K involved in the immunerelated signalling pathway was signifi cantly up-regulated (Tables S1 and 2), indicating that these immune genes might play a role in the modulation of the titre of wMel and in the immune reaction of the host to the Wolbachia. The results obtained in the present study confi rmed the reliability of the transcriptome sequencing. The RGs selected here are extensively used in other studies. Taking Actin as an example, it encodes a major component of the protein scaffold that supports the cell and determines its shape. Actin is well expressed in most types of cells and is widely used as an RG in the qPCR analysis of many different insects, including whitefl y (Zhou & Li, 2016), desert locust (Van Hiel et al., 2009), European honey bee (Scharlaken et al., 2008) and two species of Collembola (Van Hiel et al., 2009). Our results indicate that Act in and EF-1α are the most stable RGs for gDNA samples of B. tabaci, while Act in and GAPDH were the top two RGs for cDNA samples. Similarly, EF-1α is used as an RG for salmon (Olsvik et al., 2005), humans (Shen et al., 2010;Silver et al., 2006), Orthoptera (Van et al., 2009) and Hymenoptera (Horňáková et al., 2010. Nevertheless, GAPDH, the commonly used RG in many species is reported to be unstable in some studies (Tong et al., 2009;Martins et al., 2016). It is also noteworthy that RPL29 is a recommended reference gene for several conditions, including different host plants, TYLCV infection, different developmental stages and tissues . However, it is not recommended for use in studies on the artifi cial transfection of exogenous Wolbachia into B. tabaci, which further supports the notion that different reference genes should be used in different systems.

CONCLUSIONS
We developed a set of suitable RGs for studying the host insect-Wolbachia interaction using the agricultural pest B. tabaci artifi cially transfected with an exogenous strain of Wolbachia. The RGs are suitable for both gDNA and cDNA templates. A more reliable qPCR analysis of the titre of Wolbachia and the expression of functional genes will increase our understanding of the infection dynamics of Wolbachia in this pest insect, which might provide a scientifi c basis for the development of a CI-based control strategy for this pest. Fig. 2. Validation of stable RGs based on the qPCR analysis of the relative titre (gDNA) of wMel normalized by a combination of the two most stable RGs (Actin and EF-1α) and two least stable RGs (18S and GAPDH). The data are the means ± standard deviation of the three replicates. G 0 is u sed as the control in the statistical analysis using One-way ANOVA followed by SNK test at 0.05, 0.01 and 0.001 levels of probability (ns -not signifi cant; *P < 0.05; **P < 0.01; ***P < 0.001). G 0 -adults that were microinjected with the wMel strain; G 2 -2 nd generation after transfection; G 4 -4 th generation after transfection; G 6 -6 th generation after transfection; G 7 -7 th generation after transfection. Fig. 3. Validation of stable RGs by qPCR analysis of the relative expression levels of AMP/Pl3K normalized by the combination of the two most stable RGs (Actin and GAPDH) and two least stable RGs (RPL29 and ATPase). (A) AMP, antimicrobial peptide Alo-3like; (B) Pl3K, phosphatidylinositol-bisphosphate 3-kinase catalytic subunit delta isoform. The data are means ± standard deviations of three replicates. TI-G 1 is used as the control in the statistical analysis using One-way ANOVA followed by SNK test at 0.05, 0.01 and 0.001 levels of probability (ns -not signifi cant; *P < 0.05; **P < 0.01; ***P < 0.001). WT -wild-type; TI -transfected isoline; G 0 -adults emerged after microinjection with the wMel strain; G 1 -1 st generation; G 4 -4 th generation. Ribosomal protein S6 kinase beta-1 LOC109037417 1.42 a Only the genes with q value < 0.005 and |log 2 (FC)|>1 (signifi cantly regulated) are listed. FC -fold change. The same below. b The gene in bold was used in the validation of the RGs.