Selection and Validation of Reference Genes for qRT-PCR Expression Analysis of Candidate Genes Involved in Olfactory Communication in the Butterfly Bicyclus anynana

Real-time quantitative reverse transcription PCR (qRT-PCR) is a technique widely used to quantify the transcriptional expression level of candidate genes. qRT-PCR requires the selection of one or several suitable reference genes, whose expression profiles remain stable across conditions, to normalize the qRT-PCR expression profiles of candidate genes. Although several butterfly species (Lepidoptera) have become important models in molecular evolutionary ecology, so far no study aimed at identifying reference genes for accurate data normalization for any butterfly is available. The African bush brown butterfly Bicyclus anynana has drawn considerable attention owing to its suitability as a model for evolutionary ecology, and we here provide a maiden extensive study to identify suitable reference gene in this species. We monitored the expression profile of twelve reference genes: eEF-1α, FK506, UBQL40, RpS8, RpS18, HSP, GAPDH, VATPase, ACT3, TBP, eIF2 and G6PD. We tested the stability of their expression profiles in three different tissues (wings, brains, antennae), two developmental stages (pupal and adult) and two sexes (male and female), all of which were subjected to two food treatments (food stress and control feeding ad libitum). The expression stability and ranking of twelve reference genes was assessed using two algorithm-based methods, NormFinder and geNorm. Both methods identified RpS8 as the best suitable reference gene for expression data normalization. We also showed that the use of two reference genes is sufficient to effectively normalize the qRT-PCR data under varying tissues and experimental conditions that we used in B. anynana. Finally, we tested the effect of choosing reference genes with different stability on the normalization of the transcript abundance of a candidate gene involved in olfactory communication in B. anynana, the Fatty Acyl Reductase 2, and we confirmed that using an unstable reference gene can drastically alter the expression profile of the target candidate genes.


Introduction
Butterflies have fascinated mankind throughout the ages with their wing coloration, association with angiosperms and life cycle involving metamorphosis. These visually striking members of the class of insects are placed in the order Lepidoptera from which they diversified~100 million years ago and show remarkable diversity with over 17,500 known species [1][2][3]. Over the last few decades, a significant amount of work underpinning the fields of ecology, evolution, genetics, and developmental biology has been carried out with model species of butterflies. The availability of the recently sequenced genomes of Heliconius melpomene [4], Danaus plexipus [5] and Melitaea cinxia [6] has offered an extraordinary opportunity to study the expression patterns of gene families in butterflies. Besides genomic resources, tissue-specific expressed sequence tags (ESTs) libraries have also been generated for certain butterflies [7,8]. The number of butterfly genome sequences available and genome projects currently underway indicate an upward trend in the interest of researchers towards this group of organisms. Yet, an extensive study aimed at identifying suitable reference genes to normalize the expression profiles of genes of interest in any butterfly species is surprisingly still lacking.
The African bush brown butterfly Bicyclus anynana (Nymphalidae; Satyrinae) (Butler, 1879) has emerged as a model in evolutionary developmental biology since significant advances have revealed how genetic variation between individuals interacts with developmental processes to provide adaptive phenotypic variation in wing patterns, a major trait under natural and sexual selection [9][10][11]. B. anynana has also become a model in ecology and evolution, thanks to the coupling of in-depth field studies and experimental work conducted on labadapted population. This strategy has unraveled the ecological relevance of a large range of phenotypes associated, for example, with seasonal polyphenism producing adaptive phenotypes to the alternating wet and dry seasons in East Africa, life-history traits relevant to ageing research and to sexual selection [10,12]. In this regard, it has been demonstrated that B. anynana males use a male sex pheromone (MSP hereafter), produced by wing structures named androconia, which plays a key role in mate choice and sexual selection [13][14][15]. After landing behind a female, B. anynana males flicker their wings and erect their androconial hair, likely favoring the dissemination of MSP at short range over the female antennae during courtship [14]. The MSP composition has been identified and includes three chemical components [14]. Variation in abundance of the three MSP components between males informs females about male phenotypic traits associated to fitness, such as age and inbreeding status, and this in turn affects male reproductive success [15]. Interestingly, the number and position of androconia is a key taxonomic trait to distinguish between closely related species of the highly speciose Bicyclus genus (>88 currently described species) [16,17]. Moreover, the large diversity of wing chemical composition found across a subset of 32 species of Bicyclus was recently shown to have contributed to the diversification of the genus by recurrent reproductive character displacement [18].
While our understanding of the adaptive value of phenotypic variation in sex pheromone communication in B. anynana is increasing, the study of its molecular basis is still in its infancy and this holds true for olfactory communication in butterflies in general. In this regard, we decided to investigate the molecular basis of sex pheromone communication in B. anynana. For this, we needed to use quantitative real-time reverse transcription PCR (qRT-PCR) to validate the variation in level of expression of candidate genes. In this respect, reliable reference genes need to be selected and validated since non-validated reference genes can drastically alter the expression profiles of genes [19]. The present study thus aims at identifying and validating suitable reference genes in B. anynana by evaluating the mRNA expression stability profile of twelve potential reference genes under two different experimental conditions employing the widely used technique of qRT-PCR. This technique has been successfully used to identify suitable reference genes in various plants, animals and insects including for example Arabidopsis thaliana [20], Sus domesticus [21]) Drosophila melanogaster [22], Spodoptera exigua [23] Bombyx mori [24], Chilo suppressalis [25], Plutella xylostella [26], Helicoverpa armigera [27], Agrotis ipsilon [28] and Spodoptera litura [29]. As previous studies have revealed that variation in the expression profiles of reference genes are often observed in varying experimental procedures, it has become mandatory to test multiple reference genes in multiple experimental conditions to select the best reference genes for normalization studies in any given tissue set [30,31] or experimental condition [32]. Moreover, the number of reference genes needed to accurately normalize the qRT-PCR expression data must be determined on a case by case basis since using one reference gene can induce errors in the results [33]. We thus used here two algorithm-based methods, NormFinder [34] and geNorm [35] for identifying suitable reference genes, to: (a) identify the best reference gene among our 12 candidates genes; and (b) ascertain the number of reference genes needed for optimizing the normalization of gene expression levels in B. anynana. We tested our 12 candidate reference genes by sampling eight tissues of insects reared under two different experimental conditions (food stress and control feeding ad libitum). Finally, we validated the selection of our reference genes by evaluating across different tissues, the expression profile of a fatty-acyl reductase gene (FAR 2 hereafter) involved in sex pheromone biosynthesis in the Lepidoptera species B. anynana [36].

Rearing of insects
Experiments were performed with B. anynana from an outbred laboratory stock population established in the laboratory from eggs collected from a laboratory raised colony in Leiden, The Netherlands, which was itself originally established from 80 gravid females collected in Malawi in 1988 [37]. Around 200 to 400 adults per generation are reared in order to maintain high levels of heterozygosity [38]. Animals were reared in a climate room under a standard temperature regime (27°C ± 2°C), 12:12 L:D and high relative humidity (66% ± 10%) representing the wet season in the field [39]. Pupation and emergence dates were recorded, pupae were sexed and virgin females and males were kept in different cages. We applied two feeding treatments: in the control feeding treatment, all larvae and adults were fed ad libitum with maize (Zea mays) and moist banana (Musa acuminata), respectively [10]. In the food stress treatment, individuals were starved at the larval stage for 48h by rearing them in a plastic petri dish with 2 cm 3 of agar and water, while adults were fed ad libitum. We chose this larval food stress because it was shown to affect various life history traits in B. anynana including developmental time, pupal mass, fat percentage, resting metabolic rate and thorax ratio [40].

Sampling collection
We collected butterfly tissues encompassing a range of experimental conditions: three different types of tissues (wings, brains, antennae), two developmental stages (pupal and adult), two food treatments (food stress and control) and the two sexes (male and female). Eight different tissues were collected from each of the feeding treatments: eight tissue samples from the control feeding and another eight from the food stress. The tissues consisted of the pupal male wing parts containing the developing androconia (Pupae M); the pupal female entire wing tissues acting as control (Pupae F); the adult male wing parts containing the androconia (Wings M And); the adult male wing parts remaining after removal of the androconia and acting as a control (Wings M Cont); the adult female entire wings acting as a third control (Wings F); the adult male (Brain M) and female (Brain F) heads after removing the antenna, eyes and proboscis; and the pooled adult male and female antennae (Antennae M, F). We produced biological replicates of each type of tissue for both the control feeding and food stress treatments by pooling the tissues of five individuals of different ages (from 1 to 14 days-old) per sample. Bulks of different tissues collected from five individuals (three biological replicates) were used for total RNA extraction. 0.5 μg of RNA from each tissue was reverse transcribed to synthesize the cDNA.

Extraction of total RNA and cDNA synthesis
The individuals were dissected on ice directly after sacrifice and tissues were kept in RNA later (Life technologies Europe, Gent, Belgium) at -80°C according to manufacturer's instructions for total RNA extraction. Total RNA was extracted using RNeasy Mini Kit (Qiagen Benelux) following manufacturer's instructions. An optional step of DNaseI treatment was included during the RNA extraction process to degrade the residual DNA from the samples. The concentration and purity of the isolated RNA samples was measured using Nanodrop ND-1000 Spectrophotometer (Isogen Life Science) and RNA samples measuring the ratio OD 260 /OD 280 in range of 1.8-2.0 were considered for RT-qPCR experiments. The integrity of isolated RNA samples was analyzed using Aligent Bioanalyzer (Aligent Biotechnologies) (S1 and S2 Figs.). First strand cDNA synthesis was carried out for each RNA sample using Superscript First-Strand Synthesis System for RT-PCR (Life technologies Europe, Gent, Belgium) following the manufacturer's protocol in a final volume of 20 μl.

Choice of candidate reference genes and primer design
The nucleotide sequences of the ten candidate reference genes were retrieved from the transcriptome database (http://147.99.108.61:9011/ngspipelines/#!/NGSpipelines/Bicyclus% 20anynan)based on sequence homology with commonly used reference genes ( Table 1). The nucleotide sequences of the genes have been submitted to NCBI Genbank database (Accessions KM923784-KM923795). The primer sequences for two genes, eEF-1α and FK506 were retrieved from Pijpe et al [41]. For other genes, primer sequences were carefully designed using the publicly available Primer 3 software (http://bioinfo.ut.ee/primer3-0.4.0/) with the following criteria: GC content of 40-60%, Tm 60°C, length of 18-22 nucleotides and amplicon size varying between 80-120 base pairs. The secondary structures of primer sequences were checked using the software Gene Runner version 5 (www.generunner.net, version 5.0.39 Beta). Specificity of primer pairs were checked both by checking the qRT-PCR amplified products on 2% agarose gel stained in Ethidium Bromide (EtBr) and performing an in silico ePCR program [42].

Real-time quantitative reverse-transcription PCR (RT-qPCR)
RT-qPCR (or qRT-PCR) was performed in a 96-well thermocycler (StepOnePlus Real-Time PCR System, Software v2.1, Applied Biosystems-Life technologies) using QuantiTect SYBR Green PCR Kit (Qiagen, Benelux) and qPCR Master Mix plus one for SYBR Green I kit (Eurogentec) according to the manufacturer's instructions. Each reaction mixture consisted of 2.5 ng of cDNA template, 0.3 μl each of forward and reverse primers, 10 μl of SBYR green mix adjusted with nuclease free water to a final volume of 20 μl. The thermocycler program included an initial denaturation step of 10 minutes at 95°C, followed by 40 cycles of 15 second at 95°C and 1 minute at 60°C. Melt curve analysis was performed after 40 cycles by heating the samples from 65°C to 95°C to confirm the specificity of primer pairs. Amplifying primer pairs for an intron sequence on cDNA checked for genomic DNA contamination. All tested samples were technically duplicated and all experiments were performed with three biological replicates.

Data Analyses
The amplification efficiency of all primer pairs used in the study was calculated from the standard curve obtained from a five point 10-fold serial dilution series of cDNA template using the equation E (%) = (10 -1/slope -1) Ã 100 [43]. The expression levels of 12 candidate reference genes were calculated by Cycle threshold (Ct hereafter) values. The variations observed in the expression levels of reference genes were shown using the derived boxplot R-package BoxPlotR (http://boxplot.tyerslab.com).
The two most popular tools NormFinder [34] and geNorm [35] were used to evaluate the expression stability of all candidate reference genes across tissues and treatments (control and food stress). To meet the assumptions of the models used by NormFinder and geNorm, we used the method described in [44] and transformed the raw Ct values into 2 (-ΔCt) values, in which ΔCt is the difference between the Ct value of the sample tested and the minimum Ct value of the reference gene across all samples. The minimum Ct value (maximum expression level) of each gene was treated as control and was assigned a value of 1.
NormFinder uses a model based approach to assess the expression stability of different candidate reference genes by (a) measuring the overall expression level variation and (b) analyzing an intra-and intergroup variations of the reference genes. Top ranked reference genes will have minimum intra group variation and the lowest stability value. On the other hand, geNorm measures an expression stability value (M) for each reference gene by relying on the principle that the expression ratio of two ideal reference genes is identical in all samples. Thus, for every reference gene, geNorm determines the pairwise variation with all other reference genes as the standard deviation of the logarithmic transformed expression ratios, and defines the M value as the average pairwise variation of a particular reference gene with all other reference genes [35]. Moreover, Vandesompele et al [35] also indicated that more than one reference gene should be pooled in a "normalization factor (NF)" to best normalize qRT PCR expression data.
To determine the most appropriate number of reference genes to use, the pairwise variation V n/n+1 can be calculated between two sequential normalization factors, i.e. between the two most stable reference genes and as many additional normalization factors as additional reference genes were tested by stepwise inclusion of the most stable remaining reference genes. A large pairwise variation means that the added reference gene has a significant effect and should preferably be included for producing a reliable normalization factor.

Estimation of expression profile of FAR2 gene
The relative expression profile of FAR 2 gene was calculated in the eight different tissues (see materials and methods) using the 2 (-ΔCt) method. The significance of differences in gene expression levels between tissues was tested using the software R after log-transformation of the expression levels for obtaining homoscedasticity of the data across tissues, using a nested ANOVA. The independent variable, the qPCR expression levels, was compared between tissues provided that two technical replicates were produced per biological replicate and three independent biological replicates were produced per tissue type. The structure of the statistical model was thus: model <aov (log(expression level)~tissue / biological replicate / technical replicate + error (tissue / biological replicate / technical replicate)).

Results and Discussion
Selection of candidate reference genes and amplification efficiency Twelve reference genes commonly used in different organisms were selected as potential reference genes for B. anynana. The protein sequences of the reference genes from various insects were used as a query to perform a TBLASTN search in B. anynana transcriptome in order to obtain the sequences of the closest B. anynana orthologs. The selected candidate reference genes have various biological functions and include proteins involved in protein degradation (UBQL40), pentose phosphate pathway (G6PD), glycolysis (GAPDH) and ATP hydrolysis (VATPase). Also included are ribosomes encoding proteins (RpS8 and RpS18), cytoskeletal protein (ACT3), transcription factor (TBP), heat shock protein (HSP) and translation initiation factor (eIF2) ( Table 1). In addition, two reference genes previously reported in relation with ageing research in B. anynana [41], eEF1 alpha and FK506 were also included in the present study (Table 1). We paid attention to selecting genes that belonged to different functional classes, which significantly reduces the chance that genes might be co-regulated. Two primer pairs were tested for each of the 12 genes and we selected for each gene the primer pair that gave the best amplification efficiency ( Table 2). The amplification efficiency of primer pairs ranged between 90%-103% that suggests an efficient qRT-PCR system for B. anynana. The amplicon specificity of each primer pair was confirmed by the presence of a single melt curve peak (S3 Fig) and the presence of a single PCR amplicon of expected size on agarose gel electrophoresis (S3 Fig). Expression profile of the 12 candidate reference genes RT-qPCR was carried out to measure the transcript abundance of the candidate reference genes in a set of different types of tissues, developmental stages, food treatments and sexes (see Material and Methods). The potential reference genes exhibited a wide variability in expression profile under both experimental conditions (food stress and control feeding ad libitum). Overall, the observed Ct values for the 12 genes ranged from a minimum Ct value of 20.92 (eEF-1α) to maximum Ct value of 39.21 (ACT3) (Fig. 1A). The gene eIF2 showed least average transcript abundance in control and food stress treatments (Ct = 36.95± 0.26 and 36.21± 0.42 cycles, mean± standard deviation, respectively) and it showed minimum cycle of variation of 2.4 and 3.8 cycles in control and food stress treatments, respectively. The two most abundant genes were eEF-1α (Ct = 23.27± 0.24 and 24.26± 0.26 cycles, respectively) and RpS8 (Ct = 23.87± 0.42 and 24.77± 0.29 cycles, respectively) and showed moderate variation in Ct range of 3.3 and 3.1 cycles respectively (Fig. 1A). When the gene expression profiles across eight different tissues were analyzed separately for control (Fig. 1B) and food stress treatments (Fig. 1C) seven genes, namely eEF-1α, UBQL40, RpS8, RpS18, GAPDH, Act3 and G6PD showed a significant Ct variation (p value < 0.05; Ansari-Bradley Non-parametric homoscedasticity test) in food stress treatment when compared with control treatment (S4 Fig). For the remaining five reference genes (FK506, HSP, VATPase and HSP), we did not observe a statistically significant variation between the food stress and control treatments. The variation in Ct range fell between 2.4 to 6.2 cycles for all candidate reference genes for the control treatment (Fig. 1B), and increased between 3.8 to 12.8 cycles in the food stress treatment (Fig. 1C and S1 Data). Thus, no candidate reference gene had a very stable expression level across the tested tissues and treatments.

Expression stability of candidate reference genes
We used the two most commonly used statistical programs NormFinder and geNorm to monitor the expression stability of twelve candidate reference genes. Zhu et al [45] showed that these two programs are sufficient to select the appropriate reference genes. NormFinder calculated the stability values of the 12 genes and ranked the genes accordingly (Table 3). When NormFinder analysis was applied to pooled control and food stress treatments, the most stable reference gene was FK506 (M value = 0.21) followed by RpS8 and RpS18 (Table 3). In the control feeding treatment, the most stable reference gene was RpS8 (M value = 0.22) followed by FK506 and UBQL40 (Table 3). In the food stressed treatment, FK506 with an M value of 0.13 was the best choice for a reference gene, followed by RpS18 and RpS8 (Table 3). In both control and food stressed conditions, ACT3 was the least stable reference gene.
Second, we used geNorm, for which Vandesompele et al [35] strongly recommends a threshold expression stability M value of 0.5. When M values from both control and food stressed treatments were collectively evaluated, four out of twelve reference genes showed M value lower than 0.5, and RpS8 had the lowest M value of 0.35 and therefore was the best reference gene ( Fig. 2A). In the control treatment, 6 out of 12 reference genes had a M value lower than 0.5 (Fig. 2B) while in the food stressed treatment, five out of twelve tested reference genes had a M value lower than 0.5 (Fig. 2C), among which, RpS8 was each time the most stable reference gene with an acceptable M value of 0.31. Our earlier conclusion that ACT3 was the least stable reference gene was also supported bygeNorm analysis.
Thus, both NormFinder and geNorm propose the ribosomal protein RpS8 as the best reference gene under the control treatment (Table 4). Ribosomal proteins are abundant in different tissues at all different stages and therefore various ribosomal proteins have been successfully selected and validated as best reference genes in insects [23][24][25][46][47][48]. In the food stressed treatment, both software packages ranked FK506 as the best choice of reference gene. Our findings are in partial agreement with of a previous study by Pijpe et al [41] that monitored five reference genes in B. anynana: while FK506 was previously found as one of the most suitable reference gene for normalizing expression analysis data of tissues under food stress in B. anynana [41], eEF-1α that appeared to be the best choice as a reference gene in Pijpe et al [41], displayed a low ranking (fifth position) in our study (Table 4). Two reasons can explain the discrepancies observed between our   and Pijpe et al's results: first, we used a larger number of candidate genes that allowed us to find a reference gene, RpS8, that outcompeted FK506 in some conditions; second, Pijpe et al's study ranked the candidate genes using the geNorm software only, while it is strongly recommended to use both geNorm and NormFinder to identify the best reference genes [49].
Although actin proteins are used as reference genes in honeybee Apis mellifera [50] and in desert locus Schistocerca gregaria [51], our analysis ranked actin (ACT3) as the least stable reference gene. Our results confirm that actin appears to be one of the most unsuitable choices for RT-qPCR data normalization in various taxa [52]. Similarly, GAPDH that showed maximum stability as a reference gene in the moths Bombyx mori and Spodoptera exigua [25], was however not selected as a stable reference gene in our study. It appears that traditional reference genes like Actin and GAPDH may not be the best choice as normalizing genes for all organisms. Consequently, our results show that new reference genes can outperform the traditional reference genes and signifies the importance of identifying and validating suitable reference gene in B. anynana.
Although several of the 12 candidate genes showed sufficiently stable expression across different tissues to be used as reference genes in B. anynana expression analyses, we further tested whether the use of more than one reference gene could improve the normalization of the qRT-PCR expression data. We used geNorm that calculates the pairwise variation (V n /V n+1 ) between groups of reference genes that are ranked according to their expression stability across samples (see materials and methods). In both our experimental treatments (food stress and control feeding ad libitum), geNorm provided a V2/3 below the cut off value of 0.15 indicating that the use of a maximum of two best reference genes should accurately normalize the expression results (Fig. 3). For control tissues, the combination of the two reference genes RpS8 and Determination of the suitable number of reference genes for normalization. Pairwise variation (V) was calculated for (a) total tissues (b) control tissues and (c) food stressed tissues for the 12 candidate reference genes using geNorm. Dotted lines represent the recommended threshold in [35]. RpS18 (V2/3 = 0.116) could adequately normalize the expression data. Similarly, RpS18 with FK506 (V2/3 = 0.117) was sufficient for a correct quantification for tissues under food stress. Similarly, there is no increase in the optimal number of reference genes required for accurate normalization when tissues under both conditions were collectively analyzed and reference genes Rps8 with UBQL40 (V2/3 = 0.13) appeared the best for data normalization (Fig. 3). Thus, Rps8 and Rps18 were recurrently selected as part of the best pair of reference genes across our range of experimental tissues and treatments. Yet, because different pairs of reference genes were found to form the best combination in different treatments, we strongly recommend evaluating the identity of the best reference genes according to the specific experimental requirements and conditions.
Validation of the choice of reference genes using a real case study in sex pheromone biosynthesis To evaluate the effect of choosing a least or most stably expressed reference gene on the normalization of expression profile of genes, we measured the transcript abundance of the Fatty Acyl Reductase 2 gene (Ban-wFAR2 in [36]), a gene shown to be involved in the biosynthesis of the male sex pheromone component Z9-tetradecen-1-ol (MSP1 hereafter) in B. anynana [36]. We normalized the expression results for FAR2 using the best (Rps8) and least stable (ACT3) reference genes as ranked by NormFinder (Table 4). When Rps8 was used, the FAR2 transcript was significantly more expressed in male wing androconial tissues that produce sex pheromone [14,36], compared to control adult male and female wing tissues (Fig. 4A) that do not display the androconial tissues and display a significant lower abundance of MSP1 [14]. The expression of FAR2 was absent in all other tissues (Fig. 4A). Similar results were obtained when the data was normalized using the combination of the two best genes suggested by geNorm (RpS8 and RpS18; Fig. 4B). In contrast, when the least stable ACT3 candidate reference was used to normalize the expression pattern, we observed an overexpression of FAR2 in the male control tissues while FAR2 appeared less expressed in other tissues (Fig. 4C). Since androconial wing tissues are the primary sites of sex pheromone biosynthesis in B. anynana [14], the presence of very low expression of FAR2 in this tissue as compared to control tissues is unexpected. Similarly, the expression of FAR2 in the antennal tissues (Fig. 4C) is highly unlikely as there is no MSP1 in the antenna [53]. Consequently, we confirm for B. anynana that the choice of a reference gene to normalize the relative expression profile of genes of interest can alter the final outcome of the expression analysis and should be the object of a careful selection.

Conclusion
To the best of our knowledge, our work is the first in-depth study to identify and validate reference genes for qRT-PCR data normalization in a butterfly. In the present study, twelve candidate reference genes were identified and their expression profiles were measured across different tissues and conditions to identify the best reference genes in B. anynana for qRT-PCR data normalization. Ribosomal protein RpS8 was identified as the best choice of reference gene for normalizing gene expression data of genes involved in olfactory communication. Additionally, we found that actin ACT3 being the least stable should not be used as a choice for an internal control. We expect an increased interest in the lepidopteron community to assess the expression profile of gene families in B. anynana and compare it to other butterflies.