Evaluation of the expression stability of reference genes in Apis mellifera under pyrethroid treatment

Honeybees (Apis mellifera L.), which unquestionably play an economically important role in pollination and agricultural production, are at risk of decline. To study changes in gene expression in insects upon exposure to pesticides or other external stimuli, appropriate reference genes are required for data normalization. Since there is no such gene that is absolutely invariable under all experimental conditions, the aim of this study was to identify the most stable targets suitable for subsequent normalization in quantitative experiments based on real-time polymerase chain reaction in honeybee research. Here, we evaluated the expression of fifteen candidate housekeeping genes from three breeding lines of honeybees treated with pyrethroids to identify the most stable genes. The tested insects were exposed to deltamethrin or lambda-cyhalothrin, and then, changes in the accumulation of selected transcripts were assessed, followed by statistical analyses. We concluded that AmRPL32, AmACT and AmRPL13a were the commonly recorded most stable genes in honeybees treated with the selected pyrethroids.


Results
Determination of the specificity of the designed primers. In this study, we evaluated fifteen candidate genes from honeybees to check their stability in insects exposed to pyrethroid treatment ( Table 1). The main goal of the research was to identify the most stable genes that could be used as internal controls in experiments based on RT-qPCR to determine or verify differentially expressed genes in honeybees treated with pyrethroids.
First, from the GenBank database, we retrieved cDNA sequences of A. mellifera L. encoding these genes and used the data as input in Primer3 for designing the best primer pairs. By doing so, fifteen primer pairs that matched the implemented parameters were chosen, with the resulting amplicon lengths between 100 and 250 bp and the annealing temperature of the primers set at approximately 60 °C. Then, by performing a pilot experiment (end-point RT-PCR), we tested the designed primers for their specificity. The products of the RT-PCR were resolved on an agarose gel, and after staining, a single DNA band was detected for each tested primer pair (Fig. 1). No amplification products were detected in the no-template control reactions. Moreover, the results from Sanger sequencing of the cloned amplicons verified the sequence specificity of the primers used.
Next, we examined the expression rates of the tested transcripts: after each RT-qPCR, all the C T output data were grouped in a table, and a combined box plot was prepared. All of the fifteen tested genes were amplified by RT-qPCR. The highest expression rate was observed for Am18S rRNA, with a C T value of approximately 5.61. However, we omitted Am18S rRNA in further analyses because of the extremely high accumulation of rRNA in the analysed insects. The C T values of all other candidate genes ranged between 14.33 and 25.62, which were the values for AmEF1α and AmGST, respectively (Fig. 2). Importantly, the expression levels of the analysed genes were similar among the three tested breeding lines (see Supplementary Fig. S1). Next, analysis of melting curves generated during the melting stage of RT-qPCR verified the presence of a single amplification product in each reaction: the generated melting curves were sharp and symmetric (Fig. 3), indicating reaction specificity. The melting temperature of the amplicons ranged from 75.46 to 84.33 °C, which were the values for AmHMBS and Am18S, respectively. Stability analysis of candidate reference genes. In the present study, we focused on fourteen candidate HKGs of A. mellifera L. treated with two pyrethroids: deltamethrin and lambda-cyhalothrin. Our general goal was to identify the most stable reference genes in honeybees (a) regardless of breeding line and chemical compounds used for treatment and (b) individually for three breeding lines, and (c) for the two pyrethroids used for treating the insects. To achieve this goal, we used the following statistical tools: geNorm, BestKeeper, NormFinder and ΔCT. The comprehensive analysis was performed using RefFinder (https ://www.heart cure. com.au/reffi nder).
Search for the most stable reference genes for all breeding lines exposed to pyrethroids. geNorm software, installed as a Bioconductor "NormqPCR" package in R software 16 , was used as the first program to identify the stability of the predicted HKGs, and the resulting values were used to order the HKGs from the most stable (the lowest M value) to the least stable (the highest M value). The AmHMBS gene and the AmSDHA gene had the highest expression stability values (0.0597 and 0.0604 M, respectively), followed by AmTub, AmTBP, and AmRPL32. The five least stable genes were AmRPL13a, AmGST, AmRPL18S, AmGAPDH and AmEF1α (Table 2).
Next, three additional methods were used to calculate the most stable reference genes within the pool of all tested samples. NormFinder analysis classified the AmRPL32, AmHMBS, AmCHS6, AmTub and AmAct genes as the most stable genes, whereas the AmEF1α, AmGST, AmTBP, AmGAPDH and AmDORS genes were indicated as being the least stable. The ΔCT analysis showed that AmRPL32, AmHMBS, AmTub, AmCHS6 and AmAct had the highest expression stability in comparison to AmEF1α, AmGST, AmTBP, AmGAPDH and AmDORS, which exhibited unstable expression. Then, based on the standard deviation (SD) of the C T measurements, the stability values for the expression of fourteen candidate reference genes were calculated using the BestKeeper program, which showed slight differences compared to previous algorithms. The AmRPL32, AmSDHA, AmCHS6, AmHMBS, and AmRPL13a genes were identified as the most stable genes, and the AmTBP, AmEF1α, AmDORS, AmRP18S and AmGST genes were identified as the least stable genes.  (Table 4), the most stable genes were AmTBP, AmDORS, RPL32, AmHMBS and AmRP18S (according to geNorm); AmDORS, AmRP18S, AmEF1α, AMHMBS and AmRPL13a (according to NormFinder); AmDORS, AmRP18S, AmRPL32, AmGST and AmRPL13a (according to BestKeeper); and AmDORS, AmRP18S, AmEF1α, AmHMBS and AmRPL13a (according to the ΔCT method).
On the basis of the abovementioned results, the RefFinder analysis identified AmDORS, AmRP18S, AmR-PL13a, AmGAPDH, and AmEF1α as the most stable genes for the Alpejka breeding line.
Thus, the following genes were selected as the most stable genes according RefFinder calculations for the Nieska line: AmARGK, AmRPL32, AmRP18S, AmRPL13a and AmTub ( Table 5).
Analysis of the HKG stability with regard to the active substance of pyrethroid insecticide. Analysis of the influence of each insecticide used for insect treatment on HKG stability was also performed ( Table 6). All calcula-   www.nature.com/scientificreports/ tions were performed as stated above using geNorm, NormFinder, BestKeeper, ΔCT and RefFinder. In insects treated with deltamethrin, the rank order of the most stable genes was as follows: AmSDHA, AmTub, AmHMBS, AmTBP and AmRPL32 (according to the geNorm method); AmCHS6, AmRPL32, AmARGK, AmHMBS and AmRPL13a (according to NormFinder); AmRPL32, AmRPL13a, AmGST, AmAct and AmRP18S (according to the BestKeeper); AmCHS6, AmRPL32, AmTub, AmHMBS and AmARGK (according to the ΔCT method). The comprehensive analysis performed by RefFinder ranked the most stable genes in the following order for deltamethrin treatment: AmRPL32, AmCHS6, AmTub, AmARGK and AmHMBS.
The same calculations were performed to select HKGs stably expressed in A. mellifera L. exposed to lambdacyhalothrin. The C T data obtained after RT-qPCR were grouped in a table and subjected to subsequent calculations. The geNorm method indicated the following genes as being the most stable: AmTub, AmTBP, AmHMBS, AmSDHA and AmCHS6. Next, the list of the five most stable genes calculated by NormFinder included AmRPL32, AmRPL13a, AmHMBS, AmAct and AmTub. BestKeeper selected AmGST, AmRPL32, AmRPL13a, AmRP18S and AmAct as the most stable genes. AmRPL32, AmRPL13a, AmHMBS, AmAct and AmTub were identified as the most stable genes using the ΔCT method.
Determination of the minimum number of reference genes necessary for normalization. Pairwise variation analysis (V n/n+1 ) performed using the geNorm method 17 indicated that the expression of the target gene in each considered experimental variant needs to be normalized using two selected reference genes. This was indicated by pairwise variation (V) with the threshold value set at 0.15 17 . In all tested experimental variants, the V 2/3 value was lower than 0.15 (Fig. 4).
Validation of reference genes. To validate the obtained results (the indicated stable HKGs for each experimental condition), we performed an analysis of the expression of two cytochrome P450 monooxygenase (AmCYP450) genes in honeybees exposed to pyrethroid treatment. CYP450s are known to be involved in xenobiotic detoxification in insects 14 . Importantly, it was described previously that the expression of AmCY-P6AQ1 15 and AmCYP305a1 13 was influenced by insecticides. Validation experiments aimed at normalization www.nature.com/scientificreports/ of the expression of both AmCYP450 genes were performed in the following contexts: first, for the entire set of tested Carnolian honeybees exposed to pyrethroid treatment; second, for testing the effect of pyrethroid treatment on the expression of AmCYP450 in each breeding line separately; and finally, for analysing the expression of AmCYP450 separately in deltamethrin-or lambda-cyhalothrin-treated insects. As indicated earlier in the manuscript, using two reference genes is sufficient for accurate normalization of genes in pyrethroid-treated insects. For normalization of AmCYP450 expression in Carnolian honeybees, the two following HKGs were used: AmRPL32 and AmHMBS.
The results showed (Fig. 5) that expression of the AmCYP6AQ1 gene increased slightly in honeybees treated with deltamethrin (1 h and 24 h after treatment) and lambda-cyhalothrin (24 h after treatment) with a 1.35-fold change, 1.28-fold change and 1.47-fold change (all with p < 0.01), respectively. On the other hand, the expression of AmCYP305a1 in pyrethroid-treated honeybees increased over time, reaching a 4.91-fold change 24 h after treatment (p > 0.05) in the insects exposed to lambda-cyhalothrin (Fig. 5).
Next, the expression of AmCYP450 genes was validated in each breeding line individually (Fig. 6). For the Kortówka line, the AmAct and AmARGK genes were chosen; for Alpejka, AmDORS and AmRP18S were used for normalization; and for the Nieska line, the AmRPL32 and AmRPL13a genes were selected as the best normalizers for expression of the mentioned AmCYP450s.
When considering the expression of AmCYP450s in honeybees treated with pyrethroids, we observed that the expression of AmCYP6AQ1 and AmCYP305a1 in insects belonging to the Alpejka breeding line changed slightly. On the other hand, in the Kortówka breeding line, the level of expression of the AmCYP305a1 gene increased slightly 1 h and 24 h after deltamethrin treatment. For the Nieska breeding line, the expression of AmCYP6AQ1 was slightly upregulated in insects treated with both pyrethroids, whereas the expression of AmCYP305a1 was downregulated in honeybees treated with deltamethrin. These data also showed that each breeding line of tested insects responded differently to pyrethroid treatment, taking into account changes in the expression level of the AmCYP450s genes tested (Fig. 6).
Finally, the expression levels of the AmCYP6AQ1 gene and the AmCYP305a1 gene were validated separately in insects under deltamethrin treatment and lambda-cyhalothrin treatment (Fig. 7, Table 6). The active substances used in our research model have little effect on changes in the expression level of the analysed AmCYP450s genes. In particular, when analysing the effect of deltamethrin on AmCYP450 expression, a minor increase in the expression level of the AmCYP6AQ1 gene 1 h post treatment and 24 h after pyrethroid exposure (1.40-fold www.nature.com/scientificreports/ change and 1.26-fold change, respectively, p < 0.01) was indicated. Accordingly, in lambda-cyhalothrin-treated insects, AmCYP6AQ1 showed a modest, statistically significant increase in expression level 24 h after pyrethroid treatment (1.34-fold change, with p < 0.01). Similarly, the expression level of the AmCYP305a1 gene was somewhat stable over time, reaching a 1.35-fold change (with p < 0.05) in bees 1 h after treatment with lambdacyhalothrin (Fig. 7). Additionally, changes in AmCYP450s expression levels normalized against two unstable HKGs were also analysed (see Supplementary Figs. S2, S3 and S4). The use of inappropriate normalizers in differential gene expression analysis resulted in increased statistical significance at the expense of an increased error range and changes in the expression levels of target genes in individual research models (e.g., for the Nieska line, the AmCYP305a1 gene expression level 24 h after lambda-cyhalothrin treatment was almost 40 times higher than that obtained if the least stable genes were selected (see Supplementary Fig. S3). Moreover, the use of the highly unstable HKGs for validation gives different, highly discrepant results, as in the case of the Kortówka line, where after using the most stable genes, a decrease was observed in the level of AmCYP305a1 gene expression (1.19fold change), while using the least stable genes resulted in a 3.11-fold increase in the expression of a given gene (see Supplementary Fig. S3).

Discussion
To minimize both biological and experimental errors in quantitative analyses performed by means of real-time qPCR, it is important to choose the most stable reference genes for normalization of RNA input. However, this requires an individualized research approach for each analysed parameter 9,17 . One such parameter is the fitness and mortality of bees associated with commonly used insecticides, which has been extensively discussed [18][19][20][21][22] .
In this study, we investigated the expression stability of 14 candidate reference genes of A. mellifera L., belonging to Carnolian honeybees, exposed to pyrethroids. The selected subspecies of the honeybee was treated with two insecticides: deltamethrin 23 and lambda-cyhalothrin 24 . It should be remembered that honeybees of various genetic background (like the three breeding lines described in the study: Alpejka, Nieska and Kortówka) might react differently at the level of insecticide sensitivity 25,26 what can expressed at the molecular level.
Carnolian honeybees are highly adapted to nectar and climatic flow both in Poland and worldwide 27 . Analysis of all the obtained data described in this study on Carnolian honeybees under pyrethroid treatments indicated    Analysis of the expression stability of selected candidate reference genes with respect to individual breeding lines distinguished a common high-scoring gene, AmRPL13a, in terms of stability for all the tested lines. On the other hand, AmDORS, AmAct and AmTub were selected as the most stable genes in the Alpejka, Kortówka and Nieska breeding lines, respectively. Similarly common most stable genes were also observed between the Kortówka and Nieska lines (the AmARGK gene) and between the Alpejka and Nieska lines (the AmRP18S gene). The expression level of the AmARGK gene does not change after carbon dioxide narcosis in honeybee workers 28 ; however, it should be noted that the amount of ARGK protein in the antennae can vary between bee families 29 . The ribosomal genes (from the functional rRNA-coding regions) are structurally conserved and homogeneous throughout the nuclear and mitochondrial genomes in honeybees 30 and are often used as reference genes for differential expression studies [31][32][33] . In research on the effects of imidacloprid treatment on honeybees, ribosomal genes have been shown to be upregulated 13 , which means that they should be approached with caution as potential reference genes. The analyses also show the variable levels of expression of target genes relative to the AmDORS gene described in the literature 34 . Depending on the breeding line tested, the expression stability results for individual genes were classified slightly differently (Tables 3, 4, 5); therefore, in experiments, both the population and the breeding line should be determined with full accuracy to avoid statistical errors in research.
The stability ranking of HKGs in honeybees under pyrethroid treatment, when the active compounds (deltamethrin or lambda-cyhalothrin) were considered separately, showed three common most stable reference genes: AmRPL32, AmTub and AmHMBS. These results were confirmed with the data previously obtained when active substances were analysed together; however, it should again be noted that the genes were placed at different positions in the ranking order (after deltamethrin treatment: AmRPL32, AmCHS6, AmTub, AmARGK and AmHMBS; after lambda-cyhalothrin treatment: AmRPL32, AmAct, AmRPL13a, AmHMBS and AmTub). Such differences may occur due to differences in the sample sizes analysed individually. For the entire set of Carnolian honeybees, all data obtained in the experiments were taken into account. In turn, for the analysis of bees after treatment with deltamethrin or lambda-cyhalothrin, data obtained for a specific pyrethroid active substance www.nature.com/scientificreports/ treatment/exposure were taken (limiting the sample size from 108 bees up to 72 individuals). This is why the selection of the sample is such an important aspect of research related to differential gene expression 9 .
The validation of the indicated most stable reference genes showed that the selection of inappropriate normalizers, the expression of which is not stable under the conditions being tested, can significantly affect the final results of the analysis of the target gene of interest. The values may vary by up to 40 times, as was observed for the expression level of the AmCYP6AQ1 gene in the Nieska line exposed to lambda-cyhalothrin (24 h after treatment), when we compared the results obtained by using the most stable genes and least stable genes for normalization (see Supplementary Fig. S3). The statistical significance and direction of changes in the level of expression between two time points were also divergent after the selection of relatively less stable reference genes for analysis (as was the case for the Nieska breeding line, as stated earlier) (see Supplementary Figs. S2, S3, S4). Therefore, optimization of testing data by using various statistical programs is very important when studying changes in the expression of target genes relative to that of a reference gene 12 . It should also be noted that the selection of HKGs may differ if the research model assumes testing on populations, not on specific breeding lines, as presented in Tables 2, 3, 4, 5. Previous work also showed that the differences in the expression stability results for reference genes may be due to the season in which the study was conducted 35 and the stage of maturation of the tested individuals 36 . The expression levels of the AmCYP450 genes validated against the selected HKGs confirmed some behavioural observations for the developmental lines tested. Namely, slight changes in the expression of the AmCYP6AQ1 and AmCYP305a1 genes were observed for the Alpejka line (Fig. 6), in which individuals showed the highest liveliness among the three breeding lines during the experiment. Accordingly, the most considerable increase in the expression of these genes was demonstrated for the Nieska line (Fig. 6), the individuals of which tended to gather in groups and exhibited low activity (data not shown).
To summarize, regardless of the experimental conditions or tested breeding line of the examined insects, the above studies indicated the three following HKGs as reference genes to be considered, as they were classified by each analysis as being the most stable genes: AmRPL32, AmAct and AmRPL13a.
Methods insects used in the study. In this study, three breeding lines of A. mellifera L. were used, namely, Kortówka, Alpejka and Nieska, all belonging to Carniolan honeybees. The insects were taken from original hives by a beekeeper and were individually treated with a 1 µl dose of one of the following pyrethroids: deltamethrin (0.75 ml/L (4.8%)) or lambda-cyhalothrin (0.75 ml/L (4.81%)). Non-treated insects were used as a control. Then, the treated bees were gathered (6-15 insects, whereas for RNA isolation 6 insects were taken) in bee cages and collected 1 h and 24 h after treatment. During this time, the insects were kept at room temperature on a laboratory bench. Then, the insects were immediately frozen in liquid nitrogen and stored at -80 °C.