Evaluation of reference genes for gene expression analysis by real-time quantitative PCR (qPCR) in three stingless bee species (Hymenoptera: Apidae: Meliponini)

Stingless bees are generalist pollinators distributed through the pantropical region. There is growing evidence that their wild populations are experiencing substantial decline in response to habitat degradation and pesticides. Policies for conservation of endangered species will benefit from studies focusing on genetic and molecular aspects of their development and behavior. The most common method for looking at gene expression is real-time quantitative polymerase chain reaction preceded by reverse transcription (RT-qPCR) of the mRNA of interest. This method requires the identification of reliable reference genes to correctly estimate fluctuations in transcript levels. To contribute to molecular studies on stingless bees, we used Frieseomelitta varia, Melipona quadrifasciata, and Scaptotrigona bipunctata species to test the expression stability of eight reference genes (act, ef1-α, gapdh, rpl32, rps5, rps18, tbp, and tbp-af) in RT-qPCR procedures in five physiological and experimental conditions (development, sex, tissues, bacteria injection, and pesticide exposure). In general, the rpl32, rps5 and rps18 ribosomal protein genes and tpb-af gene showed the highest stability, thus being identified as suitable reference genes for the three stingless bee species and defined conditions. Our results also emphasized the need to evaluate the stability of candidate genes for any designed experimental condition and stingless bee species.

Stingless bees are generalist pollinators distributed through the pantropical region. there is growing evidence that their wild populations are experiencing substantial decline in response to habitat degradation and pesticides. Policies for conservation of endangered species will benefit from studies focusing on genetic and molecular aspects of their development and behavior. the most common method for looking at gene expression is real-time quantitative polymerase chain reaction preceded by reverse transcription (RT-qPCR) of the mRNA of interest. This method requires the identification of reliable reference genes to correctly estimate fluctuations in transcript levels. To contribute to molecular studies on stingless bees, we used Frieseomelitta varia, Melipona quadrifasciata, and Scaptotrigona bipunctata species to test the expression stability of eight reference genes (act, ef1-α, gapdh, rpl32, rps5, rps18, tbp, and tbp-af) in RT-qPCR procedures in five physiological and experimental conditions (development, sex, tissues, bacteria injection, and pesticide exposure). In general, the rpl32, rps5 and rps18 ribosomal protein genes and tpb-af gene showed the highest stability, thus being identified as suitable reference genes for the three stingless bee species and defined conditions. our results also emphasized the need to evaluate the stability of candidate genes for any designed experimental condition and stingless bee species.
A global effort to sequence genomes of all living species on Earth is in progress 1 with the aims of mitigating the impact of climate changes on biodiversity and preserving endangered species and ecosystems. In particular, the i5K consortium leads the arthropod genomics initiatives and emphasizes on sequencing insect genomes 2 . The expanded availability of genomic data during the last decade created a favorable condition for researchers to investigate a myriad of biological processes and their molecular bases, and explore the application of this knowledge to the areas of health, agriculture, industry, and ecology. In this context, more than thirteen bee species have had their genomes sequenced, among which most are from the main tribes of corbiculate Apidae: Apini (Apis cerana, A. dorsata, A. florea, A. mellifera), Euglossini (Eufriesea mexicana, Euglossa dilemma), Bombini (Bombus terrestris, B. impatiens) and Meliponini (Melipona quadrifasciata) [3][4][5][6][7][8] . Meliponini is the more diverse group of bees and comprises over 500 species of generalist pollinators 9 distributed in the pantropical region 10 . The bees of this group have a non-functional sting and thus are called stingless bees. In the last years, the wild populations of stingless bees are experiencing substantial decline that has been attributed to several factors, including habitat degradation and exposure to toxic substances such as pesticides [11][12][13][14] .

Results
Primer evaluation and expression profiles of candidate reference genes. The genomic sequences of the eight candidate reference genes of F. varia, M. quadrifasciata and S. bipunctata were aligned and primers were designed to target similar exonic regions among the three species (Table 1, Supplementary Fig. S1). The specificity of amplification of the eight candidate reference genes was initially tested for each species by conventional PCR and qPCR using a pool of cDNA samples (including all physiological and experimental conditions tested in this study) and genomic DNA. The set of primers were designed to identify possible genomic DNA contamination, as they anneal to different exons and flank intronic regions ( Supplementary Fig. S2). The fragments amplified by conventional PCR (amplicons) were visualized by electrophoresis on 2% agarose gels. Amplicons for all genes showed a single band, except for act and tbp ( Supplementary Fig. S3). Specific amplification of the genes

Stability of candidate reference genes.
To determine the stability and rank the candidate reference genes, we used geNorm, NormFinder, Bestkeeper and delta-Ct. These programs are available at the web-tool RefFinder, which also calculated a comprehensive final overall ranking based on the results of these four different algorithms. geNorm calculates expression stability value (M value) for a candidate reference gene based on the geometric mean of the SD of all studied genes in a pairwise comparison. The reference gene with the lowest M value should be the most stable gene and an M value under 1.5 is suggested by the geNorm software as a criterion for the selection of the reference gene(s) 16 . In our study, M values for all the candidate genes was lower than 1.5 for each physiological and experimental conditions and the three bee species, F. varia (Table 3), M. quadrifasciata (Table 4) and S. bipunctata (Table 5), thus indicating that all genes could act as potential reference genes. Like geNorm, NormFinder calculates stability values based on relative values, and the most stable reference genes are those exhibiting the lowest stability values 17 ; however, a cut-off value is not suggested. Bestkeeper calculates the SD value and the coefficient of variation of each candidate gene. An SD greater than 1 indicates high variation of the expression of a gene and, consequently, its instability 18 . Our results demonstrate that some candidate genes are not stable. Gene expression varied most among the tissues, organs and body parts of all the bee species that were analyzed in this study (Tables 3, 4 and 5). The comparative delta-Ct method estimates the most stable reference gene using the SD means by pairwise comparison of two reference genes. An SD below 1 indicates stable gene expression. As observed in our results obtained with Bestkeeper, not all candidate reference genes were stable, and most of the gene expression variations were verified in the tissues, organs and body parts of the bee species (Tables 3, 4 and 5). www.nature.com/scientificreports www.nature.com/scientificreports/ The comparison of the reference gene expression ranking of the four methods revealed variations in each physiological and experimental condition and bee species (Tables 3, 4 and 5). Even so, rpl32, rps5, rps18 and tbp-af were clearly among the top three stable reference genes in most analyses. It should be noted that Bestkeeper highlighted gapdh as a stable gene when we used samples of F. varia and S. bipunctata in the development condition, and also when we used F. varia males and females (sex condition) ( Tables 3 and 5). The gene ef1-α was considered stable in M. quadrifasciata and S. bipunctata treated with pesticides by NormFinder (Table 4) and Bestkeeper (Table 5). comprehensive gene expression ranking and recommended reference genes. As the four programs (geNorm, NormFinder, BestKeeper, delta-Ct) showed different ranking orders of candidate reference genes, we used RefFinder to integrate the rankings and obtain a consensus result in order to recommend the best reference gene(s) for each physiological and experimental condition, and bee species. Based on the comprehensive ranking, rpls32 and rps18 were the most stable reference genes in the various developmental stages of F. varia, M. quadrifasciata and S. bipunctata (Table 6). In females and males (sex condition), the genes rpls32 and rps18 were highly stable in M. quadrifasciata and S. bipunctata, and the genes rps5 and tbp-af were the most stable in F. varia. In the different tissues, organs and body parts of the three bee species, rps18, rps5 and tbp-af were the most stable reference genes (rps18 and rps5 in F. varia, rps18 and tbp-af in M. quadrifasciata, rps5 and tbp-af in S. bipunctata). After bacterial injection, rps18 and rpl32 were highly stable in F. varia and S. bipunctata and rps18 and tbp-af were highly stable in M. quadrifasciata. After pesticide exposure, the genes tbp-af and rps5 were highly stable in F. varia and S. bipunctata and rpl32 and rps18 were highly stable in M. quadrifasciata.
It is important to highlight that RefFinder tool does not require PCR efficiency values to calculate gene stability and some concern was raised that the outputs of this tool would be thus biased when evaluating the best reference gene(s) 42 . We used geNorm, the most affected algorithm by different PCR efficiencies 42 , to compare the results from the geNorm PLUS (qbasePLUS, version 3.0 43 ) and the geNorm from RefFinder. For all conditions and species, both geNorm PLUS and geNorm considered the same three genes as the most stable with slight differences in the position of the genes in the ranking (Supplementary Table S2-4) (an exception was tbp-af in the sexes of S. bipunctata). Thus, for the set of primers and conditions we tested, PCR efficiency is unlikely to be relevant for the identification of the most stable genes.
Influence of reference gene choice on the relative expression of a target mRNA. The four most stable and two least stable reference genes were used to normalize the expression of the immune deficiency (imd) and abaecin target genes in the bees submitted to bacterial infection. The extent of variation in imd and abaecin expression levels was then determined for each reference gene (Fig. 2). In F. varia, the relative expression of imd using the most stable genes was down-regulated after E. coli injection (tbp-af: most stable in NormFinder and rps5: most stable in BestKeeper, t-test, p < 0.05). However, when using an unstable gene, ef1-α, the result was the opposite: expression of imd was up-regulated after injection (t-test, p < 0.05). The use of the most stable  Table 3. Ranking of six candidate reference genes for Frieseomelitta varia based on their expression stability according to analyses by geNorm, NormFinder, BestKeeper, and delta-Ct during development, between sexes, in tissues, and after bacterial injection and pesticide exposure.
www.nature.com/scientificreports www.nature.com/scientificreports/ reference genes in M. quadrifasciata resulted in up-regulation of imd expression (rps5: most stable in BestKeeper, t-test, p < 0.05). In S. bipunctata, imd relative expression did not differ between control and treatment when normalized with rpl32 and rps18, which were characterized as the most stable genes by each algorithm or in the comprehensive ranking. In general, expression of imd did not change substantially after bacterial injection, but the transcriptional levels of the antimicrobial gene abaecin was highly up-regulated compared to the non-injected controls. Normalization using all selected candidate reference genes showed the same profile (Fig. 2b). These results indicate the importance of choosing stable reference genes for normalization of target genes with slight expression changes after treatments.

Discussion
The identification of gold-standard reference genes is crucial to produce reliable qPCR results. The expression of reference genes is used to correct the fluctuations in the target gene expression levels caused by technical variations in the quantity of total RNA or in the cDNA synthesis. Since there is not a universal set of reference genes that are stably expressed in all organisms and under different physiological and experimental conditions, the stability of candidate genes must be tested. Here, we evaluated the stability of six reference genes (ef1-α, gapdh, rpl32, rps5, rps18 and tbp-af) in three stingless bee species during their development, in males and females, in tissues, organs and body parts, and after bacterial injection and pesticide exposure.
Our study showed the ribosomal protein genes rpl32, rps5 and rps18 as highly stable in each of the physiological and experimental conditions for the three stingless bee species and at least one of these genes was among the most stable genes for each condition tested. The genes rpl32 and rps18 were most frequent among the genes displaying high stability and were the most stable during development (F. varia, M. quadrifasciata and S. bipunctata), between sexes (M. quadrifasciata and S. bipunctata), and after bacterial injection (F. varia and S. bipunctata) and pesticide exposure (M. quadrifasciata) ( Table 6). Ribosomal proteins are involved in translation and their genes have been found to be suitable reference genes for studies using bee species and other insects under several experimental conditions. www.nature.com/scientificreports www.nature.com/scientificreports/ The rpl32 gene is highly stable during development 19 and in nurses and foragers of A. mellifera (Hymenoptera) 21 , and in different tissues of Plutella xylostella (Lepidoptera) 44 . rps18 is among the most stable genes during the development of Tribolium castaneum (Coleoptera) 45 , Leptinotarsa decemlineata (Coleoptera) 46 and Lipaphis erysimi (Hemiptera) 47 , and in nurses and foragers of A. mellifera [19][20][21][22][23][24] , and also after bacterial infection in A. mellifera 20 and fungal infection in T. castaneum 48 . In addition to rpl32 and rps18, other ribosomal protein genes in insects are stably expressed during development, between tissues, and under different treatment conditions tested 38,44,45,49,50 .
We found the gene tbp-af as stably expressed in F. varia (between sexes and after pesticide exposure), M. quadrifasciata (among tissues and after bacterial injection), and S. bipunctata (among tissues and after pesticide exposure). The tbp-af gene codifies an evolutionarily conserved protein that composes the complex required for transcription initiation by RNA polymerase II 51 . tbp-af has also been shown to be stably expressed between tissues and after JH-treatment of A. mellifera 19 , and after exposure of Bemisia tabaci (Hemiptera) to the pesticides imidacloprid and buprofezin 52 .
Our results showed that the genes ef1-α and gapdh were the least stable genes in all tested conditions for F. varia and M. quadrifasciata. For S. bipunctata, the pair of least stable genes varied among different conditions: ef1-α and tbp-af were the least stable during development; ef1-α and gapdh were the least stable among tissues and after pesticide exposure; ef1-α and rps5 were the least stable between sexes and after bacterial injection. The protein codified by the gene ef1-α promotes chain elongation during polypeptide synthesis at the ribosome and gapdh codifies a key enzyme in the glycolytic pathway. ef1-α was also found to be the least stable gene in Bombus terrestris for target gene expression quantifications in tissues and after virus infection 26 , and also during development of the green peach aphid Myzus persicae (Hemiptera) 41 . gapdh was unstable after insecticide treatment of Plutella xylostella (Lepidoptera) 44 , during the development of Sesamia inferens (Lepidoptera) 49 , and in the green peach aphid under biotic and abiotic conditions 41 . In contrast, ef1-α was shown to be stably expressed across the tissues of A. mellifera 19 and Bombus lucorum 25,26 , and during development and between tissues of P. xylostella (Lepidoptera) 44 and during development of Diabrotica virgifera (Coleoptera) 39 . gapdh was also considered very www.nature.com/scientificreports www.nature.com/scientificreports/ stable in several studies with A. mellifera [19][20][21][22][23][24] . The different behavior of housekeeping genes among insects, and even among bee species (Table S1), underscores the need to evaluate them as reference genes for each species and physiological or experimental conditions. The high variation of Cq values observed among tissues, organs and body parts for most candidate reference genes is due to the higher Cq values obtained for ovary samples in all species (Supplementary Data S2). The reason of these variations may be due to some heterogeneity of the samples. All three species belong to a group with variable reproductive capacities in workers, ranging from workers that always activate their ovaries to sterile workers 53 . Workers of Melipona and Scaptotrigona can lay reproductive haploid eggs, which will develop into males, and throphic eggs, which are larger eggs lacking the nucleus and that serve as food source to the queen 53,54 . Frieseomelitta workers are known to lack functional ovaries and be permanently sterile, however, a recent study on the morphology of adult worker ovaries revealed the presence of normal ovarioles, suggesting that these workers may retain some capacity to activate their ovaries 53 . The variety in worker reproduction capacity in these species, as well the variety of ovary activation/degradation grade among the workers of a given species, certainly reflects in the heterogeneity of tissues in the ovaries. Heterogeneous set of samples imposes extra challenges in the search for stable reference genes 55 . Yet, our results of geNorm PLUS revealed average gene stability values (M) ≤ 0.5 ( Supplementary Fig. 10), indicating that the reference genes are stable for this set of biological samples (fat body, head and ovaries). Other candidate reference genes should be evaluated only if average M-values are >1, that indicate low reference gene stability 55 .

Condition
All candidate reference genes were tested using RT-qPCR assay to evaluate their influence on the imd gene expression in F. varia, M. quadrifasciata and S. bipunctata infected with bacteria. imd is a member of the Imd pathway, one of the main pathways of immune response in insects 56  www.nature.com/scientificreports www.nature.com/scientificreports/ case, the use of highly stable reference genes (tbp-af and rps5) evidenced down-regulation of imd expression in infected bees whereas the use of one of the least stable genes (ef1-α) pointed to up-regulation of imd expression after infection (Fig. 2). Down-regulation of imd expression was also observed in A. mellifera infected with E. coli 57 . Our experiments also highlighted that in M. quadrifasciata and S. bipunctata, imd is mildly affected by bacterial injection, thus suggesting a distinct immune response in comparison to F. varia. Different responses to bacterial infection in the recognition and signaling genes from the immune pathways have also been observed between A. mellifera and Bombus terrestris, but the up-regulation of antimicrobial peptides such as abaecin, hymenoptaecin and defensin is always observed is these bees [57][58][59][60] . Yet, up-regulation of abaecin was detected in injected F. varia and M. quadrifasciata regardless of the candidate reference gene used to normalize abaecin expression levels. This result agrees with the observations that pronounced differences in expression levels of target gene are more easily detected and less prone to be affected by the choice of reference genes 55 .
This is the first study to validate reference genes for qPCR analysis in stingless bees and provides new resources for research on this diverse and important pollinator group. Recently, more than ten bee genomes have become available 3,4,6 and allow the exploration of comparative genomic approaches to better understand the multifaceted aspects of bee biology and evolution. Yet, few studies evaluating candidate reference genes in bees were done (Table S1). Here, we used robust tools to evaluate the stability of the candidate reference genes in a comprehensive set of physiological and experimental conditions. Thus, our results have the potential to contribute to a wide range of molecular studies. In general, the ribosomal protein genes rpl32, rps18 and rps5, and tbp-af showed the highest stability in the stingless bees.

Material and Methods
Sample collection. Stingless bees. F. varia, M. quadrifasciata and S. bipunctata specimens were sampled directly from colonies located at the University of São Paulo, in Ribeirão Preto and in São Paulo, Brazil. We tested the candidate reference genes using different developmental stages, tissue (fat body), organ (ovaries) and body part (head), both sexes, and different experimental conditions, i.e., immune stimulation by bacterial injection and pesticide exposure. Details of sample collection are described below (see also Supplementary Data S2). For each species and condition, samples were stored in TRIzol ® Reagent (Invitrogen) and kept at −80 °C until RNA extraction. Whole individuals were processed for RNA extraction, except when tissues, organs and body parts were collected.
Developmental stages. F. varia, M. quadrifasciata and S. bipunctata workers were sampled throughout development including last instar larvae (defecating larvae, DL), pupae (white-eyed pupae, Pw), newly-emerged (NE) Relative expression of imd was analyzed using the six selected reference genes for qPCR normalization in bees injected with Escherichia coli versus non-injected ones (control, ct). The first and the second most stable reference genes established for at least one of the algorithms (geNorm, NormFinder, BestKeeper, and delta-Ct) are indicated by arrows. Asterisks indicate p < 0.05 (t-test). Bars represent the means and standard errors of three biological replicates. (B) Relative expression (log scale) of abaecin was analyzed using the six selected reference genes for qPCR normalization in bees infected with Escherichia coli (E. coli) versus non-infected ones (control, ct). Asterisks indicate p < 0.05 (t-test). Bars represent the means and standard errors of three biological replicates. and foragers (FOR). DL have empty intestines 61 . Newly-emerged bees (<24 h-old) were obtained from pieces of brood combs taken from the nest and kept in an incubator (28 °C) for few days. Foragers were collected in the entrance of the nest carrying pollen in the hindlegs. For each species and developmental stage, we sampled three individuals for further analysis.
Tissues, organs and body parts. Fat body, head and ovaries from F. varia FOR, M. quadrifasciata NE and S. bipunctata NE were dissected in RNase-free 0.9% NaCl solution. Each sample of fat body or head was a pool obtained from three to six individuals of F. varia or S. bipunctata. Eighteen to 22 pairs of ovaries of F. varia and 11 to 13 pairs of S. bipunctata were pooled to compose each ovary sample. For M. quadrifasciata, only one individual was used to prepare each sample due to the bigger size of this stingless bee. Pooled (F. varia and S. bipunctata) or individual (M. quadrifasciata) samples were collected in triplicates.
Female vs Male. NE females of the three stingless bee species were compared with males. NE females and M. quadrifasciata males (<24 h-old) were collected as described above for NE females. Age-specific males were obtained by marking NE males on the thorax with a paint marker, returning them to experimental hives, and collecting them after seven days (F. varia) or after eight days (S. bipunctata). Three males were collected for each species.
Bacterial injection. Ten NE workers of F. varia, M. quadrifasciata and S. bipunctata in a total of 30 bees were separated in groups of five individuals. One group of each species was injected with a bacterial suspension of Escherichia coli DH5α whereas the other group did not receive any treatment (control group). F. varia workers were injected in the thorax with 1 µL of bacterial suspension in 0.9% NaCl (5 × 10 3 bacteria/µL) using a Nanofil 10 µL syringe (World Precision Instruments, Sarasota, FL, USA) and a needle with a 0.11 mm outer diameter. M. quadrifasciata workers were injected in the abdomen (between the 5 th and 6 th abdominal segments) with 1 µL of bacterial suspension (5 × 10 4 bacteria/µL) using a Hamilton syringe and a needle with a 0.47 mm outer diameter. S. bipunctata workers were injected in the thorax with 1 µL of bacterial suspension (5 × 10 3 bacteria/µL) using a Hamilton 701 N 10 µL syringe (Hamilton ® , Reno, NV, USA) and a needle with a 0.47 mm outer diameter. The bacterial concentration used for each species was estimated based on the respective body weight. It is known that injection of 5 × 10 4 bacteria/µL is enough to activate the immune response in honey bees 62 ; thus, we used the same concentration to inoculate M. quadrifasciata, which has a body weight equivalent to A. mellifera. F. varia and S. bipunctata are about 10-fold smaller than M. quadrifasciata and the bacterial concentration used for injections was 10-fold lower. The syringe and needle types used to administer the injections were chosen according to cuticle hardness (softer in F. varia). Only injection per se is enough to trigger the honey bee immune response 60 . Thus, immune stimulation was provoked by both injection and bacterial infection. Both injected and control groups of the three species were kept in the incubator at 28 °C with food (50% sucrose solution in water) and water ad libitum. The bees were individually collected (n = 3 for each injected and control group) six hours after the bacterial injection.
Treatment with imidacloprid pesticide. To produce a primary stock solution, 10 mg of imidacloprid (Sigma Aldrich) was added to 500 µL of acetone, with a subsequent dilution in distilled water to a concentration of 1000 ng active ingredient (a.i.)/µL. This aqueous stock solution was further diluted to prepare the solutions for topical applications. NE workers of F. varia and S. bipunctata received on the thorax a topical application of 1 µL of 25.2 ng a.i./µL, a dose tested in Scaptotrigona postica 63 , whereas NE workers of M. quadrifasciata received 1 µL of 2.01 ng a.i./µL, a dose also previously tested in these bees 64 . The control bees received topically 1 µL of water on the thorax. All the bees (n = 5 for each treatment) were maintained in Petri dishes under darkness in an incubator at 28 °C with food (50% sucrose solution in water) and water ad libitum. After 24 h, the bees were collected. At this time only alive and active bees were collected (n = 3 for each treated and control group).
RnA extraction and cDnA synthesis. Extraction of total RNA was performed according to the Trizol ® manufacturer's protocol. To increase total RNA precipitation of ovary samples, 1 µg of Molecular Biology Grade Glycogen (Sigma-Aldrich) was added to the isopropanol step and samples were kept at −20 °C for 24 to 48 hours. For tissue samples (fat body, head and ovary) 1 µg of total RNA was used for cDNA synthesis. For all the other samples (DL, Pw, NE, FOR, male, NE workers injected with E. coli and their non-injected controls, NE workers treated with pesticide and their non-treated controls), we used 3 µg of total RNA. cDNA synthesis was performed using SuperScript ™ II (Invitrogen) and oligo(dT) [12][13][14][15][16][17][18] , except for the developmental stages (DL, Pw, RE, FOR) and the experiment of bacterial injection in F. varia for which SuperScript ™ III (Invitrogen) and oligo(dT) 20 were used. In both cases, cDNA synthesis was performed according to manufacturer's protocol. primer design and validation by conventional pcR. The sequences of the genes act, ef1-α, gapdh, rpl32, rps5, rps18, tbp and tbp-af were retrieved from the recently sequenced genome of F. varia (The genome assembly was submitted to NCBI under BioProject PRJNA528016) and from the genome of M. quadrifasciata deposited in the Hymenoptera Genome Database 4,65 . S. bipunctata sequences were identified in a transcriptome dataset (not published). The structural organization of open reading frames and putative exon/intron splice sites of the candidate reference genes were inferred by the annotation tool Artemis version 16.0.0 66 ( Supplementary  Fig. S2). The sequences of the genes were deposited in the NCBI data bank (accession numbers in Table 1). Intron-spanning primers were designed based on F. varia gene sequences ( Supplementary Fig. S2) using Primer3 software (bioinfo.ut.ee/primer3-0.4.0/). Melting temperatures between 60 °C and 61 °C, and amplicon length ranging from 100 to 190 bp were the restrictive parameters for primer selection. Other parameters were kept at the default setting. A maximum of two mismatches between the primer and the target sequence of M. quadrifasciata (2019) 9:17692 | https://doi.org/10.1038/s41598-019-53544-0 www.nature.com/scientificreports www.nature.com/scientificreports/ or S. bipunctata was allowed. We carefully checked the position of the mismatches in the primer sequence to avoid mismatches in the last five nucleotides of the 3′end. The only exception was the tbp-af forward primer for S. bipunctata, which contained a mismatch at the second last position 67 (Supplementary Fig. S1). All genes, accession numbers, primer sequences and amplicon sizes used in this study are listed in Table 1.
Conventional PCRs were done to check amplification from complementary DNA (cDNA) and genomic DNA (gDNA) templates. Amplifications were carried out using 2x PCR Master Mix (Promega) and 0.6 pmol of each primer, under the following PCR regime: 95 °C for 5 min; 40 cycles of 95 °C for 30 s, 60 °C for 30 s, 72 °C for 45 s; and a final extension step of 72 °C for 7 min. The PCR products were checked by electrophoresis on 2% agarose gels stained with UniSafe Dye 20.000 × (Uniscience), visualized and documented in Kodak 1D Image Analysis program, version 3.6.2 (Eastman Kodak Co., Rochester, NY).
Real-time quantitative pcR. The real-time quantitative PCR (qPCR) assays were performed using the StepOnePlus ™ Real-Time PCR System (Applied Biosystems). Amplifications were carried out in 15 μL reaction solutions containing 7.5 μL 2x qPCRBIO SyGreen Mix Separate-ROX (PCR Biosystems), 1.5 μL first-stranded cDNA (diluted 1:10), 6 ρmol/μL of each specific primer and 4.8 μL water. PCR conditions were 95 °C for 2 min followed by 40 cycles of 95 °C for 5 s and 60 °C for 25 s. The specificity of each pair of primers was checked by melting curve analysis (95 °C for 15 s, 60 °C for 1 min and a continuous raise in temperature to 95 °C at 0.3 °C/s ramp rate followed by 95 °C for 15 s). To check reproducibility, each assay was performed with technical triplicates for each of the three biological samples. PCR efficiency values (E) were calculated for each gene and bee species from the given slope after running standard curves and following the formula E = (10 (−1/slope) −1)×100. evaluation of reference gene expression stability. To determine the stability of candidate reference genes we used RefFinder 37 . RefFinder is a web-based analysis tool that integrates four algorithms: geNorm 16 , NormFinder 17 , BestKeeper 18 , and delta-Ct 68 . Based on the rankings generated by each algorithm, a weight is assigned to each gene and geometric means of the gene weights are calculated for a comprehensive final ranking. Candidate genes having lower mean weights are considered transcriptionally stable and can be used as ideal reference genes.
Validation of reference gene selection. To evaluate each candidate reference gene for qPCR assays using F. varia, M. quadrifasciata and S. bipunctata samples, we measured the expression of the imd and abaecin gene in adult workers treated with E. coli and their respective untreated controls. imd gene encodes a protein that plays a role in the Imd pathway, which controls antibacterial defense 57 . Primer sequences for imd (forward 5′AAC AAC CGA TGC AAA ACC TG 3′ and reverse 5′ TCG TTG TTT TCG GTT CAT CA 3′) and for abaecin (F. varia forward 5′-GAA GGT AAC GAC GTT TAT TTT CG-3′, reverse 5′ TGG AAA CGG ATG TCG TTG TA 3′; M. quadrifasciata forward 5′ ATG CGC GAT ATT TGC GAT A 3′, reverse 5′TTT TCG GAT TGA ATG GTC CT 3′) were designed using Primer3 software, following the same parameters described in the previous section. The relative transcripts levels of imd and abaecin were calculated for each bee species using the 2 −ΔΔCT method 69 and the six candidate reference genes. Student's t-test was performed to compare treated and control groups with significance reported for p < 0.05.