The Challenge of Stability in High-Throughput Gene Expression Analysis: Comprehensive Selection and Evaluation of Reference Genes for BALB/c Mice Spleen Samples in the Leishmania infantum Infection Model

The interaction of Leishmania with BALB/c mice induces dramatic changes in transcriptome patterns in the parasite, but also in the target organs (spleen, liver…) due to its response against infection. Real-time quantitative PCR (qPCR) is an interesting approach to analyze these changes and understand the immunological pathways that lead to protection or progression of disease. However, qPCR results need to be normalized against one or more reference genes (RG) to correct for non-specific experimental variation. The development of technical platforms for high-throughput qPCR analysis, and powerful software for analysis of qPCR data, have acknowledged the problem that some reference genes widely used due to their known or suspected “housekeeping” roles, should be avoided due to high expression variability across different tissues or experimental conditions. In this paper we evaluated the stability of 112 genes using three different algorithms: geNorm, NormFinder and RefFinder in spleen samples from BALB/c mice under different experimental conditions (control and Leishmania infantum-infected mice). Despite minor discrepancies in the stability ranking shown by the three methods, most genes show very similar performance as RG (either good or poor) across this massive data set. Our results show that some of the genes traditionally used as RG in this model (i.e. B2m, Polr2a and Tbp) are clearly outperformed by others. In particular, the combination of Il2rg + Itgb2 was identified among the best scoring candidate RG for every group of mice and every algorithm used in this experimental model. Finally, we have demonstrated that using “traditional” vs rationally-selected RG for normalization of gene expression data may lead to loss of statistical significance of gene expression changes when using large-scale platforms, and therefore misinterpretation of results. Taken together, our results highlight the need for a comprehensive, high-throughput search for the most stable reference genes in each particular experimental model.


Introduction
The term leishmaniases includes a spectrum of diseases caused by parasites belonging to genus Leishmania that affects millions of people around the world [1]. The global strategies against the disease include improving leishmaniases diagnostic tools and notification as well as providing early treatment for patients, although drug resistance is arising as a major problem. At present, there is no effective vaccine to prevent or treat leishmaniases in humans. A major limitation in vaccine development against leishmaniases is the lack of precise understanding of the particular immunological mechanisms that allow parasite survival in the vertebrate host.
The interaction of the parasite with the host cell induces dramatic changes in transcriptome patterns in the parasite, but also in the target cells [2,3]. The host withstands the infection by inducing innate immune responses, associated with changes of gene expression of numerous genes. Successful parasites can overcome the host´s immune system and colonize the tissues, and this process is deeply dependent on transcriptional reprogramming of numerous genes including those involved in basal cell processes [3]. Analysis of gene expression in infected tissues during leishmaniases will assist with the evaluation of drug and vaccine candidates, and the understanding of the immunological pathways that lead to protection and/or progression of disease.
Real time quantitative PCR (qPCR) is a highly sensitive, specific and accurate technique that allows analysis of gene expression, comparing steady-state mRNA levels between control and test samples [4]. Besides, it can be easily scaled-up in order to perform high-throughput analysis of multiple targets at the same time [5]. Many advances have been made in this field in recent years, such as the development of technical platforms for high-throughput analysis, improving mRNA quantification and internal standard selection [6,7], and the development of powerful software for analysis of qPCR data [8]. One of the drawbacks of qPCR is, however, its sensitivity, that renders the analysis susceptible to small variations in technical protocols and therefore may lead to misinterpretation of results [9,10]. In order to avoid that problem, qPCR results are generally normalized against one or more reference genes (RG) [11] to correct for non-specific experimental variation, such as differences of quantity and quality of mRNA between the samples, enzymatic activity and different transcriptional activity between tissues [12]. These RG need to be genes whose expression remains stable at various conditions across different biological groups. There are a number of studies were the selection of several traditional candidates to RG has been made in mice model using mathematical algorithms [13][14][15][16][17][18][19][20]. All these studies highlight the need to evaluate candidate RG for accurate normalization of gene expression data. The selection of RG candidates in the Leishmania-infected mice model has been traditionally focused in genes with known or suspected housekeeping roles: 18S ribosomal RNA [21,22] (Rn18s), B-actin [23,24] (Actb), Glyceraldehyde-3-phosphate dehydrogenase [25,26] (Gadph) and Hypoxanthine-guanine phosphoribosyltransferase [27,28] (Hprt). Several authors have evaluated some of them, showing that mRNA expression of genes such as B2m [29], Hprt, Pgk1 [30], Polr2a [31], Tbp [18,19] and Ubc [9] are relatively stable, although others like 18SrRNA and Actb [9] should be avoided for normalization due to high variability across different tissues or experimental conditions. The development of experiments yielding large transcriptome data sets has, jointly with qPCR studies, identified genes differing from the traditional housekeeping genes as most stably transcribed [32][33][34]. The comprehensive identification of accurate RG for normalization of gene expression data in BALB/c mice, particularly in Leishmania target organs like spleen, remains to be completed. In this study we have used three different normalization algorithms for the evaluation of 112 candidate RG genes in BALB/c mice, in order to select the most stable genes for the accurate normalization of gene expression. Our results show that, in this model and using this technical platform, some of the genes traditionally used as RG are clearly outperformed by other genes, and therefore the rational selection of RG should be strongly considered before performing large gene expression analysis.

Biological Samples
All experiments involving animals were conducted in accordance to both European (2010/63/ UE) and Spanish legislation (Law 53/2013), after approval by the Committee for Research Ethics and Animal Welfare (CEIBA) of the University of La Laguna (Permission code: CEIBA2015-0168).
Forty seven BALB/c (purchased to Charles River Laboratories, France) mice (14-15 weeks old) were used in this study. Mice were randomly separated in two groups: (i) 23 control mice and (ii) 24 mice that were infected with 10 6 stationary-phase L. infantum promastigotes via tail vein. Mice were euthanized by cervical dislocation and spleens were removed and immediately stored in RNAlater at -70°C (Sigma-Aldrich, St. Louis, USA). All efforts were made to minimize animal suffering.
Nucleic acid purity was assessed measuring OD 260/280 and OD 260/230 ratios using NanoDrop ND-1000 (Thermo). Only samples with OD 260/280 ratios between 2.10 and 2.2, and OD 260/230 ratios between 1.8-2.2 were included in this study. RNA integrity number (RIN) was determined using 2100 Bioanalyzer (Agilent Technologies, Santa Clara, United States). RIN was > 7 for all RNA samples included in this study.

Reverse Transcription and High-Throughput Real-time quantitative PCR (RT-qPCR)
Reverse transcription was performed with random hexamers using the High Capacity cDNA Reverse Transcription Kit (Life Technologies, Carlsbad, CA) according to manufacturer's instructions. Briefly, 10 μL of total RNA (200 ng/μL) were added to the 2X reverse transcription mix for a total of a 20-μL reaction. Thermal cycler conditions were set to run for 2 h at 37°C, and for 10 min at 75°C followed by a rapid cooling to 4°C. Real-time PCR was performed using QuantStudio 12K Flex Real-Time PCR System (Life Technologies, Carlsbad, CA). Custom TaqMan OpenArray Real-Time PCR Plates were custom-designed, and included 112 Gene Expression Assays organized in 48 subarrays. cDNA samples were loaded into the plates using an OpenArray AccuFill instrument (Applied Biosystems) according to the manufacturer's protocols. Each sub-array was loaded with 2.5 μl 2X GeneAmp Fast PCR Master Mix (Applied Biosystems), 1 μl 5X TaqMan OpenArray Remix (Applied Biosystems), 0.3 μl RNase-free water, and 1.2 μl cDNA. The thermal cycling protocol was to manufacturer defaults (95°C for 15 seconds, 60°C for 1 minute, for 40 cycles). All reactions were performed in triplicate. Quant-Studio 12K Flex Real-Time PCR System (Life Technologies, Carlsbad, CA) calculates Cq values using an algorithm that takes into account the efficiency of each individual curve, called C rt method [35]. The C rt method sets a threshold for each curve individually that is based on the shape of the amplification curve, regardless of the height or variability of the curve in its early baseline fluorescence. The method first estimates a curve that models the reaction efficiency from the amplification curve. It then uses this curve to determine the relative threshold cycle (C rt ) from the amplification curve, that eliminates the need of "conventional" Real-Time PCR for calculating the efficiency of each reaction. Therefore, Cq values produced by this platform are already corrected for the efficiency of the amplification.

Assay Design and Selection of Targets
One hundred and twelve genes were included in this study, all of them available as TaqMan Gene Expression Assays (Life Technologies). In order to evaluate their performance, six candidate RG were selected after a literature search of the most commonly used RG for gene expression assays in BALB/c mice: Hprt, Ubc, B2m, Pkg1, Polr2a and Tbp. 106 extra genes were included in the array, all of them involved in different mechanisms of immune system: cytokines, chemokines and their receptors, transcription factors, prostaglandins, costimulatory molecules, enzymes and others. The complete list of genes is shown in S1 Table.
geNorm. geNorm normalization was performed in qBase PLUS software (Biogazelle, Gent, Belgium). The geNorm algorithm evaluates RG candidates using two quality parameters: the stability value (M) and the coefficient of variation of the normalized reference gene expression levels. M value is the average pairwise variation of a particular gene with all other control genes included in the analysis. geNorm allows eliminating the worst-scoring candidate RG and recalculating M values for the rest, considering those with lowest M values (M < 0.5) the most stably expressed across the data set. geNorm also determines the number of RG that are required for optimal normalization using V parameter, calculating the pairwise variation (V n/n+1 ) between the normalization factors of the samples. Generally, V < 0.15 is considered the threshold value below which an additional reference gene is not required for accurate normalization [36]. In brief, the Cq value of every sample, by triplicate, was loaded in qBase PLUS software, with set cut-off values of M < 0.5 and V < 0.15.
NormFinder. NormFinder [37] is a Microsoft Excel-based application available online (http://moma.dk/normfinder-software) that identifies the optimal normalization genes among a list of candidates. This tool ranks the candidates genes depending on their expression stability, considering the variations in expression levels between sample subgroups (intergroup variation) as well as within the subgroups (intragroup variation) assigning a stability value for each candidate gene. [38] is a website tool available online (http://fulxie.0fees.us/?type= reference) which integrates the four most common normalization algorithms: NormFinder, BestKeeper [39], geNorm and ΔCt [40] method. This tool ranks every candidate RG based on the geometric mean of the four methods mentioned above.

Normalization of Gene Expression
qBase Plus software was employed to obtain Relative Quantity (RQ) and Normalized Relative Quantity (NRQ) values from the whole data set, following manufacturer's instructions [41]. Normalization was performed using two different pairs of RG: Il2rg+Itgb2 and Polr2a+Tbp.

Evaluation of Reference Genes and Statistical Analyses
Gene expression data (NQR) was compared for three different genes (Cxcl10, Ifng and Tnf). Statistical analyses were performed using the non-parametric Mann-Whitney test for pairwise comparison between control and infected groups, using SPSS software, version 20 (IBM). Uncertainty in data accuracy was calculated as described by Remans et al., 2014 [10].

Results
Spleens from a collection of 47 BALB/c mice, divided in control mice (n = 23) and L. infantum infected-mice (n = 24) were analyzed by RT-qPCR in order to determine gene expression of 112 different genes. After determination of the quality and integrity of RNA and cDNA, and in order to validate its performance in RT-qPCR, it was verified that none of our 47 template cDNAs failed to show amplification with all 112 Taqman Assays, indicating that all of them have enough quality to produce at least one RT-PCR result. In fact, the average of positive RT-PCR amplifications was 97 Taqman Assays/sample (out of 112 possible), very high, taking into account QuantStudio 12K Flex Real-Time platform is a semi-automated system performing thousands of reactions at once. Normalization of gene expression in this massive data set (15792 data points: 112 qPCR reactions performed in triplicate using RNA from 47 mice) was performed using geNorm, NormFinder and RefFinder algorithms independently. As a first step, 41 genes had to be excluded from the analysis. The three algorithms used are sensitive to missing data, meaning that if a single Taqman assay did not yield amplification in one mouse, that gene must be excluded for RG identification across the whole data set. Only 1516 qPCR reactions (out of 15792, less than 10%) did not amplify, but that meant that 41 genes could not be included in the analysis to identify the most stable reference genes. Therefore only 71 genes were accepted for the normalization process. The selection of the most stable RG was performed initially using only control (non-infected) mice, in order to evaluate the performance of the candidate RG (the full ranking is included in S2 Table). In a second stage, gene expression stability in Leishmania-infected mice group was analyzed (the full ranking is included in S3 Table), as well as using the whole set of mice (the full ranking is included in S4 Table), in order to test the robustness of the evaluation.
Expression stability of reference gene candidates in control BALB/c mice Cq data was loaded into qBase PLUS software and geNorm was implemented. According to geNorm algorithm, only Hprt, from the typical RG group, showed an M-value below the threshold (M< 0.5), unlike Ubc, B2m, Polr2a, Pkg1 and Tbp, whose M values were higher than 0.5 (Fig 1A). This result indicates that all these typical RG (Ubc, B2m, Polr2a, Pkg1 and Tbp) do not meet the criteria to be used as reference genes in gene expression assays. Moreover, 24 genes related to different immune mechanisms showed better stability than classical RG, reflected as M-values  Bold letters indicate those genes that were identical among the top-ten by the three methods. Italics indicate "traditional" reference genes according to literature.
doi:10.1371/journal.pone.0163219.t001 below 0.5 (Fig 1A). geNorm algorithm ranked Il18bp and Il10rb as the most stable genes (Table 1), enough for optimal normalization according to V parameter (Fig 1B). NormFinder analysis of this set of mice ranked Myd88 and Il2rg as the most stably expressed genes (Table 1). Again, the ranking of stability values for the collection of classical RG (Fig 1C), showed that only Hprt was ranked among the top-20 genes (7/20), and two of them (Polr2a and Tbp) were among the 10-worst scoring genes (S2 Table). This same analysis performed with RefFinder algorithm ranked Hprt and Stat4 as the most stable candidate reference genes (Table 1). Only Hprt and Ubc were ranked among the top 20 genes, in first and thirteenth position respectively (Fig 1D).
Although the genes selected for normalization for the group of control mice by the three algorithms are different, five of the top 10 genes (Itgb2, Il2rg, Il10ra, Hprt and Myd88) are identical regardless of the selection algorithm ( Table 1), indicating that combinations of these genes are probably the best choice for normalization of gene expression data in experiments using BALB/c mice model.

Expression Stability of Reference Gene Candidates in BALB/c Mice Infected with Leishmania infantum
In order to meet the criteria that mark good reference genes, their expression must remain constant at various biological conditions across different biological groups. For that reason, the same analysis was performed on gene expression data obtained from a group of Leishmaniainfected mice, and also from the whole set of mice, both control and infected with L. infantum.
geNorm analysis showed that only two commonly used reference genes, Hprt and Pkg1, show M < 0.5, and therefore are acceptable as RG, ranked 8 th and 15 th respectively; in contrast, expression stability of Ubc, B2m, Polr2a and Tbp are all above the threshold, hence their use is not recommended (Fig 2A). geNorm ranked Il6st and Itgb2 as the most stably expressed genes ( Table 2), enough for optimal normalization according to V parameter (Fig 2B). NormFinder analysis also agreed to only rank Hprt and Pkg1 (among the candidate RG) in the top 20 genes, 7 th and 14 th respectively (Fig 2C), far from the genes whose expression is most stable in the Leishmania-infected mice group, Itgb2 and Il2rg (Table 2). In contrast, Il10rb and Tgfbr1 were the top-ranked genes for normalization of gene expression according to RefFinder algorithm ( Fig 2D). The commonly used reference gene Pkg1 only ranked 10 th in this ranking.
The gene expression of this set of genes is apparently more homogeneous in this group of mice, since 7 of the genes ranked in the top-10 are identical for the three algorithms: Il2rg, Itgb2, Stat6, Il10rb, Tgfbr1, Il6st and Tgfb1 (Table 2). Interestingly, Il2rg and Itgb2 stand out as two robust reference genes for gene expression analysis in spleens of BALB/c mice, regardless of infection with an intracellular parasite whose target organ is spleen, among others.
The stability of Il2rg and Itgb2 as preferred reference genes in this experimental situation was challenged using a heterogeneous collection of mice, i.e. both healthy control and Leishmania-infected mice. The disparity between these two groups must be reflected in the stability values of the genes. Any gene that shows good stability values in "control" group, "infected" and "control+infected" mice will provide a very good reference for accurate normalization of gene expression in this large data set. geNorm analysis of the whole set of mice revealed that 19 genes showed M < 0.5, only one of them being a traditional RG: Hprt, ranked 7/19 ( Fig 3A). According to this ranking, Itgb2 and Stat6 are the most stable genes (Table 3), and enough for an optimal normalization ( Fig  3B). It must be pointed out that Il2rg ranks 4 th in this analysis. NormFinder analysis revealed Il2rg and Il6ra as the genes with highest expression stability in this ranking ( Table 3). The best scoring classical RG is Hprt, as low as 23 rd (Fig 3C). RefFinder algorithm was also used to determine the best genes for normalization in the whole data set (Fig 3D). Again, Hprt is the only  Bold letters indicate those genes that were identical among the top-ten by the three methods. Italics indicate "traditional" reference genes according to literature.
doi:10.1371/journal.pone.0163219.t002  commonly used RG that ranks among the top-20 (6/20). In this ranking, Itgb2, Stat6 and Il2rg were most stable genes (Table 3). Four genes, Il2rg, Itgb2, Stat6 and Il10ra, ranked among the top-10 by the tree algorithms used in the whole set of mice. It is remarkable that Il2rg and Itgb2 are ranked among the best scoring candidate reference genes for every group of mice and every algorithm used, which is indicative of their robustness as reference genes in gene expression analysis in this experimental model, despite its intrinsic heterogeneity.

Evaluation of Il2rg and Itgb2 as Reference Genes
In order to evaluate the pair of reference genes selected using the three algorithms, Il2rg and Itgb2, these genes were used to normalize the whole data set, and then compared to the results after normalization with a pair of RG genes chosen upon their extensive use in the literature, Tbp [18,19] and Polr2a [31]. Hence, normalized gene expression data (NRQ) was compared among control and infected groups for three different genes (Cxcl10, Ifng and Tnf) whose expression is widely analyzed in Leishmania-infected BALB/c mice, and are well-known to be overexpressed in this model. These genes code for cytokines that play crucial roles during infection and are overexpressed after Leishmania infection [42][43][44][45].
Normalization with either set of RG revealed an increased gene expression for Cxcl10, Ifng and Tnf in Leishmania-infected mice, as expected upon literature. However, normalization using Il2rg and Itgb2 revealed statistically significant gene expression differences among infected vs control animals for all the genes tested (Cxcl10 (p < 0.01), Ifng and Tnf (p < 0.5) after Mann-Whitney tests), whereas normalization with Tbp and Polr2a did not (Fig 4). This finding emphasizes the relevance of choosing the right set of RG in gene expression experiments in order to produce reliable results.

Discussion
Identifying the most stable RG for normalization of gene expression analysis under varying experimental conditions has always been a challenge for scientists. During the last few years is becoming increasingly evident that, in order to obtain reliable gene expression data, the adequate selection of RG is key. This problem has been generally addressed by using one or more genes with known or suspected housekeeping roles, but numerous studies are reporting that transcription of these genes can fluctuate considerably under varying experimental conditions and tissues [9,46]. Moreover, most of these papers only test 5-15 candidate RG genes to identify those with the highest stability, overlooking others, not generally described as RG that could outperform the classical RG. In this paper we evaluated the stability of 112 genes using three different algorithms: geNorm, NormFinder and RefFinder in spleen samples from BALB/ c mice under different experimental conditions (control and Leishmania-infected mice). Despite the discrepancies in the stability ranking shown by the three methods, due to differences in their mathematical approach, most genes show very similar performance as RG (either good or poor) across this massive data set (>15000 qPCR reactions).
Six classical RG were initially selected as candidate RG upon their extensive use as RG for qPCR experiments in the literature: B2m [29], Hprt [47,48], Pgk1 [30], Polr2a [31], Tbp [18,19] and Ubc [9]. Surprisingly, classical RG like B2m, Tbp and Polr2a were constantly ranked in our experiments among the less stably expressed genes under these experimental conditions; therefore, we do not recommend the further employment of these genes as normalization controls in RT-qPCR analysis in this model, neither for control nor for Leishmania-infected BALB/c mice.
Only Hprt was ranked among the 10 most stable genes for most of our analyses, both in control and infected mice, and consequently is the only classical RG that we would consider using in future experiments. This finding is in agreement with several papers [49][50][51][52] that evaluated classical RG candidates in mouse intestine, adipocytes, hepatic steatosis and brain respectively. On the contrary, in other studies in mice muscle [53] and non-B lymphocytes [9] the use of Hprt gene as RG was not recommended due to considerably variation of gene expression. These discrepancies only highlight the need for a comprehensive screening of adequate RG for each experimental design in order to obtain reliable results.
Our search for the best RG in spleen samples of BALB/c mice revealed that Il18bp + Il10rb was the best combination of genes according to geNorm, Myd88 + Il2rg for NormFinder and Hprt + Stat4 for RefFinder. Despite this apparent disagreement, five of 10 the top-ranked genes are identical regardless of the algorithm: Itgb2, Il2rg, Il10ra, Hprt and Myd88. The agreement between the algorithms is even higher when analyzing Leishmania-infected BALB/c mice: seven out of the 10 the top-ranked genes are the same: Il2rg, Itgb2, Stat6, Il10rb, Tgfbr1, Il6st and Tgfb1; being Il6st+Itgb2 the best combination for geNorm, Itgb2+Il2rg for NormFinder and Il10rb+Tgfbr1 for RefFinder.
Given the interest in transcriptome changes during Leishmania infection in spleens of BALB/c mice, finally both groups (control and infected) were taken together to identify the most stably expressed genes. geNorm ranked Itgb2+Stat6 as the best combination of RG, NormFinder chose Il2rg+Il6ra and RefFinder proposed Itgb2+Stat6. Four of the best 10 RG proposed by the three algorithms are the same: Itgb2, Stat6, Il2rg e Il10ra.
It is remarkable that, despite the natural heterogeneity of this model, Il2rg and Itgb2 are always ranked among the best scoring candidate reference genes for every group of mice (control/Leishmania-infected/altogether) and every algorithm used, which is indicative of their robustness as reference genes in gene expression analysis using a high-throughput platform in this experimental model. It may seem surprising, given that most of these genes are related to immune response and have never been described as RG in the literature. Il2rg encodes a transmembrane protein that is a common subunit of several interleukin receptor complexes. Itgb2 encodes an integrin beta chain, which combines with multiple different alpha chains to form different integrin heterodimers, cell-surface proteins that participate in cell adhesion as well as cell-surface mediated signaling. Although not "traditional" housekeeping genes, both have roles related to homeostasis of immune system. After the analysis of this extensive collection of data, it is clear that these genes present higher expression stability in this model than the commonly used RG, and should be considered a suitable alternative. A similar finding was described by Müller et al., (2015) [54] in a very different model, pathogen-infected tomato leaves. After a genome-wide microarray analysis and RT-PCR, they identified a set of very stably expressed genes that were never described for that purpose before, and were always evaluated as more stable than the traditional housekeeping genes included for comparison.
Finally, in order to demonstrate the effects of using "traditional" vs rationally-selected RG for normalization of gene expression data in our model, changes in expression of Cxcl10, Ifng and Tnf genes in control/Leishmania-infected mice was determined after normalization using the pair of RG identified by our screening, Il2rg + Itgb2, or two classical RG, Tbp +Polr2a. Cxcl10, Ifng and Tnf gene expression is widely described to increase during leishmaniases in BALB/c mice [42][43][44][45]. After initial Leishmania infection of macrophages and dendritic cells, these cells produce proinflamatory cytokines and chemokines, such as TNF-α and CXCL-10, that in turn recruit new host cells (natural killer cells (NK), monocytes and neutrophils) to the infection site [55,56]. Similarly, NK cells produce IFN-γ in an attempt to activate macrophages to phagocityze parasites. Normalization with Tbp+Polr2a failed to reveal statistically significant gene expression differences among infected vs control animals for all the genes tested. In contrast, normalization using Il2rg + Itgb2 showed differential gene expression for Cxcl10, Ifng and Tnf as expected.
Taken together, our results highlight the need for a comprehensive search of the most stable reference genes in each particular experimental model. The use of "traditional" or non-evaluated reference genes in gene expression experiments may drive scientists to extract flawed conclusions from their research, which is particularly relevant when using high-throughput platforms like QuantStudio 12K Flex Real-Time PCR System. In our model of Leishmania infected BALB/c mice, extensively used for preliminary testing of vaccine candidates against leishmaniases (reviewed in de Oliveira et al., 2009 and references therein) [57], classical RG like B2m, Tbp and Polr2a have shown very poor stability as RGs in this large-scale gene-expression platform, and therefore we cannot recommend their use as normalization controls. In addition, our results demonstrate the advantage of a comprehensive, high-throughput approach to identify suitable reference genes for transcriptomic studies. In this particular experimental model (spleen samples of Leishmania infected BALB/c mice) our approach has identified and evaluated Il2rg + Itgb2 as the best combination of RG, but it is unclear whether this results can be extrapolated to other samples or animal models. Therefore, scientists planning large qPCR gene-expression experiments should make all efforts to identify the best RGs for their particular experimental conditions using a comprehensive high-throughput approach. Besides, validation of the selected RGs in each and every experimental approach is another fundamental issue that should be addressed in order to obtain reliable results from these highthroughput experimental designs.
The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [58] and are accessible through GEO Series accession number GSE80709 (https:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE80709).
Supporting Information S1