Assessment of Aedes albopictus reference genes for quantitative PCR at different stages of development

Members of the Aedes genus of mosquitoes are widely recognized as vectors of viral diseases. Ae.albopictus is its most invasive species, and are known to carry viruses such as Dengue, Chikugunya and Zika. Its emerging importance puts Ae.albopictus on the forefront of genetic interaction and evolution studies. However, a panel of suitable reference genes specific for this insect is as of now undescribed. Nine reference genes, namely ACT, eEF1-γ, eIF2α, PP2A, RPL32, RPS17, PGK1, ILK and STK were evaluated. Expression patterns of the candidate reference genes were observed in a total of seventeen sample types, separated by stage of development and age. Gene stability was inferred from obtained quantification data through three widely cited evaluation algorithms i.e. BestKeeper, geNorm, and NormFinder. No single gene showed a satisfactory degree of stability throughout all developmental stages. Therefore, we propose combinations of PGK and ILK for early embryos; RPL32 and RPS17 for late embryos, all four larval instars, and pupae samples; eEF1-γ with STK for adult males; eEF1-γ with RPS17 for non-blood fed females; and eEF1-γ with eIF2α for both blood-fed females and cell culture. The results from this study should be able to provide a more informed selection of normalizing genes during qPCR in Ae.albopictus.


Introduction
Aedes albopictus (Ae. albopictus) is the most widely-travelled member of the Aedes genus [1]. Due to increased trade and ease of transcontinental travel [2,3], they are now found abundantly in not only the tropical Asian countries from which they originate, but also temperate Asia, Europe, the Americas, Australia, and Africa [4]. Though less notorious than Aedes aegypti (Ae. Aegypti) as a vector of arboviruses, Ae. albopictus is acknowledged as an efficient vector of at least 22 viral strains [5]. It also carries a multitude of cross-species infecting bacteria [6,7]. Throughout the course of history, this mosquito has been deemed responsible for a number of major outbreaks of diseases normally related to Ae. aegypti, such as Dengue [8,9] and Chikugunya [10]. Ae. albopictus is additionally susceptible to Zika, with high dissemination and transmission rates [11]. Nonetheless, it is typically presumed that where Ae. aegypti a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 populations are rampant, Ae. albopictus would assume a secondary role to the former where the propagation of a virus is concerned. The geographical overlap and subsequent competitiveness between the two species is in fact often considered a strategy for Ae. aegypti control. However, Ae. albopictus appears to be more acclimatisable, and in more recent times has displaced Ae. aegypti as the dominant species in various locales [12]. As reports of a particular strain of Chikugunya evolving to be more transmissible by Ae. albopictus have emerged [13], worries regarding the possibility of other arboviruses also exhibiting such plasticity with regards to transmissibility by Ae. albopictus has grown and thus, understanding the host-virus interactions within the species as well as its efficiency as a propagator are areas of study gaining much interest.
An influx of gene expression studies is thus anticipated in the near future. Since it was first introduced in the 1990s, quantitative real-time polymerase chain reaction (qPCR) has become the main technique of choice for such purposes. This method requires minimal nucleic acid quantities compared to more its more traditional RNA-quantification compatriots, in addition to bearing the added advantages of being faster and more reproducible [14]. Its vast popularity has also made ready a variety of commercially accessible reagents and automated platforms. However, inconsistencies in protocol, template quality, as well as enzymatic efficiencies brought forth by such diverse product-availability have contributed to irregularities in qPCR data interpretation [15,16]. Normalization of data against a reference gene is therefore critical. This involves comparing the ratios of expression levels of the target gene against that of the selected reference gene(s) [17]. The importance of such an exercise is twofold; normalization not only compensates for variations in starting cDNA quantities amongst samples, but also lends the significance in differences seen from quantification data a higher degree of confidence [18,19]. Reference genes such as TUBB and GAPDH are examples of genes utilized heavily for this purpose. The commonality of this practice is borne out of the assumption that a reference gene is expressed consistently throughout the cell cycle. However, this is no longer considered true. No one gene has been found to be stably expressed under all developmental and experimental conditions [20][21][22][23]. A reference gene may in fact be specific towards a protocol and/or biological sample. Algorithms such as BestKeeper, geNorm, and NormFinder have since been developed to identify the best-fit reference gene to use with consideration to one's unique experimental variables [17,24].
In Ae. albopictus, ACT [25] or RPS7 [26] are often applied as the normalizing reference gene. However, the application of a solitary gene for this task has been dismissed by many as a flawed approach. The gene validation algorithms as aforementioned oftentimes recommend that at the very least, a combination of two candidates should be used together to achieve proper normalization [27][28][29]. This is especially crucial for organisms subjected to a wide array of environmental changes throughout its development. Like Ae. aegypti, Ae. albopictus is a holometabolous organism undergoing complete metamorphosis, and spends half its lifecycle as an aquatic organism. Minute diet and temperature differences are known to affect growth and development patterns of larvae and adult Ae. aegypti [30][31][32]. Expression levels of certain genes in female adults of both Ae. aegypti and Ae. albopictus have also been shown to depend on whether the individual has partaken a blood-meal [33,34]. As a number of these genes do appear to affect the mother-offspring transmission of viruses as well as the ability of the mosquito to carry viruses as adults, there is thus a need for the stability of candidate reference genes to be validated at different points of development to ensure robustness in gene expression normalization for this pathologically important species.
Here, we comprehensively evaluated nine candidate reference genes i.e. RPS17, RPL32, eEF1-γ, eIF2α, PGK1, PP2A, ILK, ACT and STK at fourteen points of development as well as C6/36 cell culture. Adults are further divided by sex and blood meal status. The three algorithms of BestKeeper, geNorm, and NormFinder were used to analyse their stability and to rank the candidate genes in order of favourability for usage within any biological sample. The suitable reference gene(s) suggested through the outcomes of this study can be applied for normalization of qPCR data for whole organism Ae. albopictus tissue at multiple developmental stages as well as cell culture.

Rearing and sample collection
Dried viable eggs of VCRU-lab strain Ae. albopictus were obtained from Vector Control Research Unit, Science University of Malaysia (USM). Roughly 500 were submerged in dechlorinated water at any one time. First instar larvae (1L) samples were collected immediately upon hatching. Hatchlings were maintained in relative humidity and natural light conditions at 28˚C in plastic containers throughout the entire span of their aquatic lifecycle. Rearing water is changed every other day. Feed of approximately 1.0g mixture of crushed dog biscuits and pulverized chicken liver was administered daily. Second (2L), third (3L), and fourth (4L) instar individuals were collected successively. Pupal samples were a mixture of an equal number of individuals at first, second, and third day of pupation. Newly-eclosed adults were transferred to cages. They were maintained on 10% sucrose solution. Combined adult samples collected comprised of equal numbers of males and nonblood fed females aged 1 to 10 DAE. Adults of the same age range were also gender-separated in sampling. Food was removed from 5 to 7 DAE adults a full day before blood-feeding on lab-strain mice. Mosquitoes were returned to normal breeding conditions for females to lay eggs. Embryos were collected at three-hour intervals up to 12 hours; six-hour intervals from 12 to 24 hours; and 24-hour intervals from 24 to 72 hours after they were laid. C6/36 cells cultured in Gibco 1 L-15 media supplemented with 10% FBS, 1% Pen-Strep and 10% Tryptose Phosphate Broth (ThermoFisher Scientific, USA) were harvested at maximum confluency.

RNA extraction and quality assurance
This study adheres to the Minimum Information for Publication of Quantitative Real-Time PCR guidelines or MIQE. The amount of tissue collected per bioreplicate are as follows: 500 eggs per embryonic sample; 50 individuals per first or second instar larval sample; 35 individuals per third larval sample; 20 individuals per fourth instar larvae or pupal sample; 20 individuals per adult sample; and 2ml of 4 th day culture per cell sample. All samples were immediately stored -20˚C in TRIzol 1 reagent (Invitrogen™, Ambion™, Life Technologies). Total RNA extraction was done within five days of collection with a protocol previously described for mosquito tissue samples. Culturing media was removed from c6/36 samples prior to RNA extraction as described by Abcam 1 . Extracts were quantified on the Quawell 1 Q3000 UV Spectrophotometer (Quawell Technology, Inc., California). The acceptable A260: A280 value was set between 1.75 and 2.05. All extracts used showed clear 18S banding and minimal smearing in 1.0% agarose gel, and were kept at -20˚C for the duration of the experiment.
were designed on the Primer3 software (bioinfo.ut.ee/primer3-0.4.0/). The Kalign nucleic acid alignment software (http://www.ebi.ac.uk/Tools/msa/kalign/) was used to specify for regions spanning exon-exon boundaries. Restrictive parameters for primer selection were: melting temperatures between 59.0˚C and 61.0˚C, GC content between 40 and 60%, nucleotide length between 18 and 24, and amplicon length of between 150 to 225 bases. Other settings were kept at default. Singular amplicon product and mRNA specificity was confirmed in silico on the Sequence Manipulation Suite website (http://www.bioinformatics.org/sms2/). All genes, accession numbers, primer sequences and amplicon size used for this study is listed in Table 1.

Reverse transcription and qPCR
Reverse transcription was carried out with 1.5μg total RNA in 30μl reactions using the iScript Reverse Transcription Supermix (Bio-Rad Laboratories, California; cat. no. 1708840) according to manufacturer's protocol. Pooled undiluted cDNA from all eleven developmental stages were serially diluted to the factor of 5 (1:1, 1:5, 1:25, 1: 125, 1:625) for standard curve generation. All qPCR runs were performed on the BioRad CFX96 qPCR platform. Optimum qPCR reactions were in 10μl reactions of iTaq™ Universal SYBR 1 Green Supermix (Bio-Rad Laboratories, California; cat. no. 1725120),~10ng total cDNA, and 500nM each of forward and reverse primers. The standard protocol is initial denaturation at 95˚C for 2.30mins, followed by 40 cycles of denaturation at 95˚C for 15s, annealing at 59˚C for 15s and extension at 72˚C for 20s. A melting-curve analysis with a temperature range of between 65˚C and 95˚C immediately followed amplification. All samples were quantified in technical triplicates. Expression levels were recorded as cycle threshold (Ct). Efficiency values (E) were calculated according to the equation: E = (10 [-1/slope] -1) X 100. were utilized for selection of best candidate gene. The BestKeeper algorithm is an Excel program which generates a ranking through repeated pairwise correlation and regression analysis of a gene against all the other tested candidates. Raw data of Ct values (annotated as CP) and PCR efficiency of the primers were used to determine the correlation between each candidate gene and the index, expressed in the form of a coefficient of determination. For geNorm and NormFinder, raw data was converted into linear values relative to the lowest Ct recorded for each candidate gene. In geNorm, the stability of a gene is assessed through the consistency of its expression ratio across all samples. The software generates both a stability value i.e. M, and a pairwise variation value i.e. V. M represents the average variation in transcript levels of a gene in comparison to all other candidate genes, achieved through a repeated process of stepwise exclusion commencing from the least stable gene. Pairwise variation estimates the effect of including another gene [17] sequentially as per the established M-value rankings through the formula of V n /V n +1. A threshold of 0.15 is set; a V value below this would mean that an additional reference gene would not improve normalization. NormFinder is a mixed-effects model statistical analysis which estimates the stability value of a gene as a function of the approximate expression variation it would impose onto the target gene data during normalization [37]. The lower this value is, the less variation one would introduce to a normalization exercise should the candidate gene be used as a reference. It also estimates the variation between sample subgroups of the sample set. The BestKeeper vs. Pearson correlation coefficient value, geNorm M value, and NormFinder stability value are perceived as weightage. Geometric means i.e. central tendencies of these weightages for a candidate gene forms the basis for generation of a consensus ranking.

Primer pair evaluation of candidate reference genes
The expression patterns of an upwards of sixteen housekeeping genes were initially observed from previously reported RNAseq data [38,39]. Seven genes including Ribosomal Protein L34 (RPL34), α-tubulin, β-tubulin, RNA Polymerase II (RNAPII), 18S, TATA-Box Binding Protein (TBP), and Ribosomal Protein S7 (RPS7) were eventually excluded due to any one or more of the following factors: (a) lack of introns or exon-exon boundaries enabling mRNA-specificity; (b) unsuitably low or high Ct values, which compromises sensitivity; and (c) poor primer design (see S1 File). Nine shortlisted candidates progressed through to the next stages of analysis. They can be divided into four functional classes: (i) ribosomal and protein-production genes: RPS17, RPL32, eEF1-γ, and eIF2α; (ii) metabolism-related gene: PP2A; (iii) signal-transduction genes: PGK1, ILK and STK; and (iv) structural integrity gene: ACT. A standard curve for each primer pair was generated with pooled cDNA serially diluted by a factor of 5. Primer quality is based on the efficiency (E) and linear regression coefficient (R 2 ) values as observed from the amplification of cDNA at very high to very low concentrations. All recorded acceptable E values between 95.0 and 117.1%. R 2 values range from 0.979 to 0.998 (Table 1). Amplification specificity was displayed through the production of a singular peak in melt-curve analysis, and confirmed on a 2% agarose gel (S1 Fig). PCR products were also sequenced and confirmed to have an alignment of at least 95% to the predicted gene region. All sequences are obtainable from GenBank's BankIt depository with the accession numbers KY199533 to KY199541.

BestKeeper analysis
BestKeeper estimates the standard deviation (SD) value of each candidate gene from raw Ct numbers. An SD>1 signifies that the variations in expression of a gene within a sample of Assessment of Aedes albopictus reference genes for quantitative PCR at different stages of development the same origin are high, and thus indicating its instability. Our data demonstrated that not all candidates were stable across all samples (see S1 Table). ILK and PGK1 frequently appeared to be unreliable. Both gave an SD above the acceptable threshold in 48 to 72 hour embryos, larval samples from second to the fourth instar stages, pupal samples, and bloodfed females. PGK1 were additionally unstable in 24 to 48 hour embryos and cell culture. STK displayed instability in embryos 24 to 72 hours in age as well as second instar larvae; PP2A in 48 to 72 hour embryos, second instar larvae, and blood-fed females; eEF1-γ in in third instar larvae and pupae; ACT in third instar larvae and cell culture; and eIF2α in third and fourth instar larval samples, pupae, and blood-fed females. Both ribosomal-linked genes i.e. RPS17 and RPL32 were stable in all developmental time points. The SD of a gene factors into its position on a ranking based on the value given as BestKeeper vs Pearson correlation of coefficient. The closer this value is to 1, the greater the reliability of the gene. Here, those carrying SD values of above 1 are relegated to the bottom of the table regardless of its correlation of coefficient value. BestKeeper most recommends ACT for 24 to 48 hour embryos, second instar larvae and blood-fed adult females; eEF1-γ for 18 to 24 hour embryos; eIF2α for each of the 18 to 24 hour and 48 to 72 hour embryonic stages along with cell culture samples; ILK for 0 to 3 hour, 3 to 6 hour, and 9 to 12 hour embryos; PGK1 for 6 to 9 hour embryos; PP2A from the fourth instar larval stage through male adulthood; RPL32 in third instar larvae; and finally RPS17 during the first instar larval period as well as non-blood fed females. Rankings are shown in Table 2.
Genes highlighted red had an SD of over the recommended stability indicator of 1. BestKeeper values (not shown here; please refer S1 geNorm analysis geNorm provides two assessment outcomes. The first is an average expression stability score as symbolized by M. The higher the M-value of a gene, the less stable it is perceived to be. This value should fall below 1.5. None of the candidate genes exceeded this in any developmental stage ( Table 3). As considerations for a pairing of genes for normalization purposes is incorporated into the algorithm, two are viewed as equally stable or applicable for any single stage. An ILK and PP2A pairing is suggested for 0 to 3 hour embryos; RPL32 and RPS17 for the 3 to 6 hour and 9 to 12 hour embryonic stages as well as blood-fed females; PGK1 each with RPS17 for 6 to 9 hour embryos, with PP2A for the 12 th to 18 th hour, and with ACT for 18 to 24 hour embryos along with the first instar larval stage. RPS17 is to be paired with ILK for 24 to 48 hour embryonic samples; with eEF1-γ for fourth instar larvae and both sexes in non-blood fed adults; and with STK for pupal samples. Second instar larval samples is best normalized by an ACT and eEF1-γ pairing. For each of the remaining developmental time points i.e. 48 to 72 hour embryos, third instar larvae, and cell culture, RPL32 provides fairest equilibration when paired with eIF2α, PP2A, and eEF1-γ, respectively. Rankings are summarized in charts as provided by the software (Fig 2). The second outcome from geNorm estimates the effect of a gene addition event [17]. The proposed cut-off value, denoted as V, is 0.15. Our data shows that all samples of larval and pupal tissue, as well as those of 48 to 72 hour embryos and non-blood fed female adults, should require more than two candidate genes to be properly normalized (Fig 3).

NormFinder analysis
Similar to geNorm, the data utilized by this program is based on relative values. The algorithm produces a stability value for each gene when compared to others within the group. A lower value indicates greater stability. The program does not make suggestions for a cut-off value [37]. Rankings generated are summarized in Table 4. eEF1-γ is deemed the best reference gene for the 6 to 9 hour embryonic stage and both male and blood-fed female adults; PGK1 for nonblood fed females, 9 to 12 hour and 12 to 18 hour embryos, and first instar larval samples; PP2A for 0 to 3 hour and 12 to 18 hour embryonic tissue as well as pupal samples; RPL32 for 48 to 72 hour embryos, cell culture, and larval samples of the second through fourth instar; RPS17 for 24 to 48 hour embryos, and finally STK for those aged between 3 to 6 hours. Despite the frequency at which PGK1 is top-ranked, it is also the most commonly disrecommended candidate (n = 4/17). In the absence of group identifiers, it is presumed that the two genes with the lowest stability value within a sample set would provide the best combination for two-reference gene normalization strategies [40].

Consensus list of reference genes
Consensus rankings are obtained through geometrically averaging the weights assigned to a candidate gene (in the form of stability values from geNorm and NormFinder, and a function of 1-((BestKeeper vs. Pearson correlation coefficient value) from BestKeeper) by the three  programs. All genes are included regardless of BestKeeper SD values. Results are summarized in Table 5. The three top-ranked genes in the consensus list for any developmental stage are typically considered to be most reliable [37]. However, as these combinations vary greatly from one developmental stage to the next, reliability across sample types is also assessed on the basis of overall frequency at which a gene appears amongst the top-three. eEF1-γ and RPL32 are the most reliable, with a frequency of 0.157 and 0.137, respectively.

Discussion
The mosquito is intensively studied for a number of reasons: the central role of several of its species as propagators of disease-causing organisms notwithstanding, the insect's wide-travelled nature and historically paralleled evolution alongside humans have also provided insights into adaptive organogenesis and development [41,42]. Ae. albopictus is especially interesting to study as it is both a potent carrier of arboviruses and the most invasive species of mosquito to have emerged in recent years [4,43]. As recognized pests, much of the research surrounding  [44][45][46][47]. Regardless, comparative assays and subsequent investigations into how a single pesticide could affect these two genetically similar organisms vastly differently could be the key to understanding the genetics of mechanisms of resistance, and eventually offer a solution for overcoming the growing problem of pesticide resistance in Ae. aegypti. These studies more often than not require qPCR. The method is hugely popular as it is robust, powerful, and fast. As it is widely utilized, many considerations for obtainment of assured data is also locked in place [28]. When these guidelines are followed, qPCR data is often viewed as highly reliable. Such a consideration is the normalization of expression levels through the application of reference genes [48,49]. These are typically chosen from a group classified as 'reference genes'. However, current opinion has shied away from the assumption that only such genes need apply as references. In fact, it is highly likely that genes most suitable  are often utilized, usually singularly and with little thought for the potential instability of these genes in respect to experimental variables. These practices are unadvisable; not only is the utilization of a single gene for normalization no longer accepted [50], but a misinformed selection of reference genes are known to lead to false positives or false negatives [15]. It is thus of utmost importance that this lack of information is promptly addressed. The challenge lies in the unannotated state of the Ae. albopictus transcriptome, which reduces confidence in the identity of genes of interest within the species. Nonetheless, through high-confidence sequence alignments in respect to the thoroughly annotated A. aegypti transcriptome, a group of candidate genes were chosen for validation within the context of this study. Ultimately, we hope that the findings here could act as a guide towards more careful selection of reference genes in qPCR assays involving Ae. albopictus. A total of fifteen pre-determined stages of development in Ae. albopictus were sampled. A subgroup of female adults at 24-hours post-blood meal, a point where RNA production is doubled [51], as well as C6/36 cells are also included. Stability of nine candidate genes within each developmental stage is separately evaluated by the BestKeeper, geNorm, and Normfinder tools. The three gave largely differing results, though a level of congruence was seen amongst top three rankings in a majority of the stages considered. The gene most frequently observed within this group was eEF1-γ. It appeared eight times out of a possible fifty-one. The gene outperformed RPL32 (n = 7/51), as well as PP2A, ACT and PGK1 (each n = 6/51). In contrast, STK was most frequently found in the bottom-three rankings (n = 8/51), followed closely by ILK and RPL32 (each n = 7/51). Given that RPL32 was equally prominent in top three and bottom three rankings, and as inconsistencies in how most of the candidates were ranked across the consensus board are seen, it is more advisable to categorize performance on a tissue-type basis. As we are also wary of preferences in the scientific community to limit reference genes to two, the following suggestions will be composed of the pair with the best showing within each tissue-type grouping.
For embryo-derived samples up to 24h post-oviposition, PGK1 and ILK is the most recommendable normalizer pairing. Both encode for kinases. Along with STK, they are chosen as candidates as previous validation studies have found members of this category of genes to be reliable references [52,53]. Kinases act as signal transducers. As information relay occurs around the clock, they are in constant demand, and this translates into the constitutive expression of their genes. However, the 'reusable' nature of these enzymatic proteins limits their numbers within the cell at any one time. This is hypothesized to assert a degree of expressional stability and subsequently, the applicability of this group of genes in normalization. Regardless, as most reference genes are regulated, their levels were shown to be affected by biological needs. It is thus important to note that both PGK1 and ILK performed relatively poorly outside of the 0 to 24 hour embryo tissue-type.
Greatest reliability amongst candidate genes in embryos aged between 24 and 72 hours shifts towards ribosomal genes. This trend carries through into the larval and pupal stages as well. RPL32 is prominent within the top-three rankings throughout these developmental periods. We therefore tentatively suggest the pairing of this gene along with RPS17 for all tissues along this span of time. Ribosomes are the protein-generating machines of the cell system. Its two subunits, the 40S and 60S ribosomal proteins, are built from a total of over 80 components. Whereas RPS17 is a component of the 40S subunit, RPL32 is one of 60S'. Both are common internal control genes and have displayed stability when utilized in cell culture [54][55][56], plants [57,58] and mammalian tissue [59,60]. In insects, RPL32 [61] and RPS17 [62] were shown to be reliable references regardless of tissue origin during bodily development. Cellular activity of the evolving mosquito speeds up after the 24 th hour within embryos. Our raw data shows expression of either gene to be high (Mean Ct value; RPL32 = 16.77, RPS17 = 16.61). When only the larval stages are considered, these values are further reduced, indicating even greater expression levels. This is the point of development whereby growth is most robust. Consequently, demand for proteins to facilitate the process is at an all-time high. Nonetheless, the kinetics of ribosomal activity is such that their components can dissociate and re-associate as needs dictate. The balancing act between availability, demand, and need here may altogether incur a stabilizing effect upon ribosomal genes and thus, their inherent usability as controls for these tissue types. This coincidentally explains why RPL32 and RPS17 are reasonably stable within rapidly growing cell culture as well. As C6/36 cells are also of larval origin, we are doubly confident with suggesting these two genes for normalization in context of Ae. albopictus cell culture.
Within adult tissues, another participant of protein synthesis emerged as the best reference gene. Eukaryotic Elongation Factor 1 Gamma or eEF1-γ is a subunit of the protein eEF1, which delivers aminoacylated-tRNAs to the ribosome during translation. However, the overall contribution of this highly versatile protein within the cell is much more widespread, potentially holding roles in nuclear export, proteolysis, and apoptosis, amongst others [63]. Our evaluation here marks eEF1-γ as superior to eIF2α as a candidate reference gene. At the outset of this study, our expectations were such that these two would be comparable as like eEF1-γ, eIF2α is also a multifunctional player of protein synthesis [64]. It is unclear to us as to why the latter is the more volatile gene. Nonetheless, these two should provide proper normalization when utilized for expression within cell culture and blood-fed adult female tissues. For adult males and females, eEF1-γ is best paired with STK and RPS17, respectively. We had also hypothesized an association between eEF1-γ and ACT, as eEF1 is known to effect actin polymerization and subsequent cytoskeletal formation [65,66]. True enough, the two were found closely placed in the consensus rankings in eight out of the seventeen sample subsets here, though unexpectedly this occurred mostly under circumstances where neither showed reliability as a normalizer.
An additional form of result supplied by geNorm is a pairwise variation (V) evaluation ( Fig  3). The addition of a suitable reference gene during target gene normalization is expected to reduce the normalization factors (NF) value. However, geNorm dictates that if the stepwise inclusion event is valued at below 0.15, subsequent addition events will not positively affect normalization outcomes. The V values for our study show that in many cases, the combination of two genes is sufficient. Interestingly, the sample subsets which appear to benefit from stepwise inclusion all belong to the 'rapid-growth' group i.e. 48 to 72 hour embryos through to pupae. Where the inclusion of none of the genes evaluated would satisfy the 0.15 threshold value, it was proposed that three would be the ideal number to use during qPCR [67,68]. For this reason, in addition to the pair of genes suggested for usage on a tissue-type basis above, we recommend the addition of a third gene most top-ranked for an individual subset as seen Table 5. Nonetheless, this observation also alludes to there being a panel of reference genes more suitable for this period of development than the ones evaluated here.
In Aegypti, we were able to recommend a two-gene combination of RPS17 and ACT as a normalization strategy across all tissue-types on the basis of their inherent stability throughout development [62]. For Albopictus, eEF1-γ and RPL32 instead are the genes most commonly seen in top-three rankings within the consensus. Despite this, eEF1-γ is also seen as the least reliable gene in two stages i.e. 9 to 12h embryos and pupae. Similarly, RPL32 is as equally 'stable' as it is 'unstable'. Therefore, we cannot assuredly advise on a two-gene combination which would be suitable enough for general use regardless of tissue-type and experimental conditions. The overall middling performance of ACT is another surprising outcome of the analysis, given the gene's history as a ubiquitously used control in not only qPCR, but also proteomics. Its stability within the developing Ae. albopictus was however only on par to the equally middling PP2A. We had included the phosphatase in this study given its verified usability as a reference gene during development in plants [69][70][71]. In animal tissue, it is less commonly studied for this purpose. Our results suggest that although PP2A could be applied effectively as a reference for six tissue types, the lack of obvious trends to its stability renders the gene unpredictable. This emphasizes the need for a re-evaluation of suitability of genes as references on a species-by-species basis. Improvements to be made in the future include broadening the classes of genes evaluated as well as putting variation-contributing factors such as viral infection or environmental stresses under parallel consideration. In vastly researched species such as Drosophila, a large panel of reference genes has also been shown to benefit normalization as sample size and experimental complexity grows [72][73][74].

Conclusion
This is the first study of its kind in Ae. albopictus. Through the algorithms of BestKeeper, geNorm, and NormFinder, we recommend the implementation of two-gene combinations to provide satisfactory normalization for the described developmental stages and C6/36 cell culture on the basis of tissue-type. Based on consensus rankings, the proposed combinations are PGK and ILK for early embryos (0 to 24 hours post-oviposition); RPL32 and RPS17 for late embryos (24 to 72 hour post-oviposition) as well as all four larval instars and pupae samples; eEF1-γ with STK for adult males; eEF1-γ with RPS17 for non-blood fed females; and eEF1-γ with eIF2α for both blood-fed females and cell culture. Application of an additional gene during the span of development from 48 to 72 hours post-oviposition embryos to the pupal stage comes highly recommended. These findings will benefit normalization practices in Ae. albopictus, and may additionally serve as a resource for screening reference genes in closely-related insects.
Supporting information S1 Fig. PCR products in 2% agarose gel. In flanking lanes are 100bp ladder. (PDF) S1 File. Expression level comparisons between candidate genes, and discarded primers. Table 1 is a summary of candidate gene expression levels and performance from previous publications; Table 2 is a list of primers eventually excluded from the study. (DOCX) S1