{'fraction_non_alphanumeric': 0.03720424997454089, 'fraction_numerical': 0.024440748158457516, 'mean_word_length': 4.941912061315047, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gene Expression Variation of Candidate Endogenous Control Genes Across Latitudinal Populations of the Bank Vole (Clethrionomys glareolus)

The bank vole has become an important species for the study of gene expression changes underlying evolutionary adaptation, within-host dynamics and resistance to pathogens, or response to pollutants. RT–qPCR is an optimized method for the rapid and accurate measuring gene expression, but it relies on the use of reference genes (RGs) as endogenous controls to normalize mRNA levels of the target genes. Contrary to the common practice, the expression level of the housekeeping genes traditionally used as RGs cannot be assumed stable across different species and experimental settings. Furthermore, compared to model laboratory species, there are potentially additional sources of variation when collecting gene expression data from natural populations with unknown genetic and environmental backgrounds. Thus, rigorous determination of RGs in natural populations of the bank vole was necessary to facilitate gene expression studies in this emerging model species. Therefore, we evaluated the expression variation of 10 RG candidates in spleen samples of bank voles spanning a broad latitudinal range across Europe, using four approaches. Ube2d2a, Ppia, and Tbp were consistently identified as the least variable genes, followed by Rn18s, Ywhaz, and Rplp0. The combinations Ube2d2a and Ppia or Ube2d2a and Tbp were identified to be sufficient for the normalization. In contrast, the traditional housekeeping genes Actb, Gapdh, Sdha, and Tuba1a displayed large expression variation and are not recommended as internal controls. Our results underscore the need of a systematic evaluation of appropriate RGs to each particular experimental system and provide a useful starting point for further studies.


INTRODUCTION
Regulatory variation in gene expression is increasingly being recognized as a major determinant of phenotypic variation (Lai et al., 2008;Fraser, 2013;Sandkam et al., 2015;Lin et al., 2017). However, the role of regulatory changes in the evolution of natural populations has not been firmly established, in part because the majority of gene expression studies have been conducted with model organisms in controlled laboratory settings (Nuzhdin et al., 2004;Fraser et al., 2010;Lin et al., 2017). Such studies explore the effect of genotype and physiological state of the organism on gene expression, but their generality is limited by the restricted range of genotypes and phenotypes available in laboratory populations (Lai et al., 2008). There is therefore great, but theretofore underutilized potential for understanding evolution and other biological phenomena such physiology, immunology, parasitology and toxicology, by the study of gene expression in natural populations (Townsend et al., 2003).
The real time quantitative reverse transcription-PCR (RT-qPCR) is an optimized method for the rapid and accurate measuring of gene expression, which provides means to compare the expression levels of target genes in individual samples, including validation of high-throughput gene expression profiles (VanGuilder et al., 2008). To obtain reliable patterns of gene expression in a given setting using RT-qPCR, it is important to minimize the influence of technical or experimentally induced variation (Pfaffl, 2004;Huggett et al., 2005). Therefore, normalization of target gene expression is typically performed by using one or multiple reference genes (RGs) as an endogenous control (Bustin, 2002). Most often, genes involved in maintaining cellular structure and homeostasis, the so-called housekeeping genes (Eisenberg and Levanon, 2013), serve that purpose, because they are expected to be expressed at constant levels across tissues (Pfaffl, 2004). However, the use of such genes without prior evaluation can lead to inaccurate results and erroneous conclusions about real biological effects (Bustin, 2002;Kozera and Rapacz, 2013). Further, compared to laboratory animals, there are potentially additional sources of variation in natural populations (Huang and Agrawal, 2016), due to genetic background differences or environmental heterogeneity (Harrison et al., 2012;Lenz, 2015), and careful selection of RGs is therefore critical.
The bank vole (Clethrionomys glareolus) is a small rodent with broad distribution in Europe (Deffontaine et al., 2005;Filipi et al., 2015), ranging approximately from 40 • N in Calabria in Italy to 69 • N in Lapland in Finland and Norway, between where it spans a broad variation in altitude and climate (Johannesen and Mauritzen, 1999;Torre and Arrizabalaga, 2008). The bank vole displays one of the most complex geographical patterns of population structure in any small mammal studied to date (Pedreschi et al., 2019), likely as the result of isolation and adaptation of populations during and after the glacial periods of the Pleistocene (Deffontaine et al., 2005;Kotlík et al., 2014;Filipi et al., 2015). While suitable glacial refugia for the bank vole have persisted in southern Mediterranean regions of Europe and in the vicinity of the major mountain ranges (such as the Carpathians), facilitating population divergence and accumulation of genetic diversity, recurrent founder events during end-glacial colonization of northern areas resulted in a reassortment and/or reduction of diversity away from the refugial areas . The broad range and complex history made the bank vole one of the key mammals used in genetic studies of the response to climate change at the end of the last glaciation (Wójcik et al., 2010;Colangelo et al., 2012;Kotlík et al., 2018;Marková et al., 2020) including local adaptation to environmental heterogeneity (Kotlík et al., 2014;Strážnická et al., 2018). Additionally, bank vole is an important zoonotic reservoir of various emerging pathogens such as Puumala orthohantavirus (PUUV) (Khalil et al., 2019;Madriéres et al., 2019) Ljungan virus (Hauffe et al., 2010) and cowpox virus (Franke et al., 2017), with significant differences in the level of host immune-related gene expression reported between bank voles from PUUV endemic and non-endemic regions of Europe (Guivier et al., 2014).
Over the past years, several traditional housekeeping (e.g., Actb, Gapdh, Rn18s, or Tuba1a) have been used for RT-qPCR data normalization in the bank vole in order to assess gene expression patterns implicated in the susceptibility to PUUV infection (Dubois et al., 2018;Müller et al., 2019), in the cellular response to the chronic exposure to ionizing radiation (Mustonen et al., 2018) or environmental heavy metal contamination (Świergosz-Kowalewska and Holewa, 2007;Zarzycka et al., 2016;Mikowska et al., 2018), and to measure PUUV viral RNA load in bank vole tissues (Guivier et al., 2014). However, the common RGs may show high expression variability in different tissues or experimental settings and other genes may thus present suitable alternatives (Hernandez-Santana et al., 2016). In a recent study (Dvořáková et al., 2020), two RGs (Ppia and Rn18s) showing low expression variance across several individuals were used to reveal the possible role in environmental adaptation for the gene expression regulation of bank vole globin genes (Dvořáková et al., 2020). However, to date, there has been no report on a rigorous statistical identification of stably expressed RGs in the bank vole or a related species. Therefore, we selected a set of commonly used as well as non-traditional RG candidates (Supplementary Table S1) and evaluated their expression variation among bank vole samples chosen to represent the species' geographical range and variation in population history across Europe (Figure 1), by using a suite of mathematical and statistical algorithms, i.e., geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), BestKeeper , and the comparative C q method (Silver et al., 2006). The results support that traditional RGs may be outperformed by other candidates. Our study thus underscores the importance of evaluating RGs for each experimental system and provides a useful starting point for further studies in the bank vole.

Samples
The study used spleen samples (n = 21) of male (n = 10) and female (n = 11) bank voles from 13 localities across the species' broad latitudinal range in Europe (Figure 1 and Supplementary Table S2). The samples were available from previous genetic Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 1 | Geographical origin of the bank vole samples (n = 21). Map shows the geographic origin of the samples in the three latitudinal zones: southern refugial (n = 6), which includes the areas of glacial refugia in the Mediterranean region, northern apex (n = 7), which includes Great Britain and Fennoscandia as the apex areas of the northward end-glacial colonization; and central (n = 8), which has a mixed history of glacial survival and post-glacial colonization. Geographically proximate sampling sites may overlap (see Supplementary Table S2 for details). studies and represented the three latitudinal zones identified based on the combined information from phylogenetic analyses of mitochondrial DNA sequences (Filipi et al., 2015) and clustering of genomic data : southern refugial (n = 6), which includes the areas of glacial refugia in the Mediterranean region, northern apex (n = 7), which includes Great Britain and Fennoscandia as the apex areas of the northward end-glacial colonization; and central (n = 8), including samples from the area between southern and northern regions, which has a mixed history of glacial survival and post-glacial colonization (Filipi et al., 2015;Marková et al., 2020).

RNA Isolation and Quality Control
Spleens stored in RNAlater (Qiagen, Hilden, Germany) at −20 • C were carefully examined and only tissue of medium size with no apparent necrotic foci or cysts were used. Approximately 0.05 mg of tissue was homogenized on ice with a T10 ULTRA-TURAX homogenizer (IKA, Staufen, Germany). Total RNA was extracted using the RNeasy Protect Mini Kit (Qiagen), followed by RNase-Free DNase I treatment. The RNA concentration was measured with a NanoDrop TM 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts, United States) and the RNA integrity evaluated with an Agilent 2100 Bioanalyzer using the Agilent RNA 6000 Pico Kit (Agilent Technologies, Santa Clara, CA). The resulting RNA of each samples was diluted to the final concentration of 5 ng/µl, aliquoted and stored at -80 • C until the analyses.
All reactions were carried out in triplicates and non-template controls were included in all RT-qPCR runs. Amplifications were performed using 0.1 ml strip tubes and caps (Qiagen) placed in the 72-well rotor. The RT-qPCR program consisted of reverse transcription at 50 • C for 30 min, initial denaturation at 95 • C for 15 min, followed by PCR cycles consisting of denaturation at 94 • C for 15 s, annealing at temperature specific for each set of primers (Supplementary Table S3) for 15 s and extension at 72 • C for 20 s; and final extension at 72 • C for 5 min. Fluorescence data were acquired at the end of each extension step. The specificities of the RT-qPCR products were verified by a melting curve analysis tool of the Rotor-Gene Q Series Software, version 2.3.1, and by gel electrophoresis on a 1.5% agarose gel and stained with GoodView Nucleic Acid Stain (SBS Genetech, Beijing, China). The amplification of the desired PCR product was confirmed by Sanger sequencing (Macrogen, Amsterdam, Netherlands) and BLAST comparison against the National Center for Biotechnology Information (NCBI) 1 database.

Data Processing
The fluorescence data were processed with the Rotor-Gene Q Series Software v. 6.1 (Corbett Research), which uses the second derivative of raw fluorescence values (the slope of the amplification curve) to help calculate both the efficiency (E) of each reaction and the take-off point at which the exponential phase of amplification begins, taken to represent the quantification cycle (C q ) value (Mccurdy et al., 2008;Bustin et al., 2009;Camacho Londoño and Philipp, 2016). For comparison, quantification cycle values were also determined by manually setting a threshold line in the linear fluorescence scale, hereafter referred to as C t values to differentiate them from the C q values obtained by the second derivative method (Mccurdy et al., 2008). As RefFinder (Xie et al., 2012) only takes quantification cycle values without the possibility to input E values, corrected C q and C t values were calculated by using the formula C q × (logE/log2). The C q and C t values were transformed to linear scale relative expression quantities (RQ) required by geNorm (Vandesompele et al., 2002) and NormFinder (Andersen et al., 2004), by using the formula E (control Cq−sample Cq) (Livak and Schmittgen, 2001;Mccurdy et al., 2008), with the sample with the lowest mean C q or C t value used as a control (De Spiegelaere et al., 2015). Relative quantity calculations based on the C t data were also performed by assuming 100% efficiency for all individual reactions by setting E at 2 (Livak and Schmittgen, 2001). Outliers, that is replicates whose amplification efficiency or relative expression clearly deviated from the sample mean, were not included in the calculations.
The C q and C t values were used with BestKeeper  while the RQs were an input to geNorm (Vandesompele et al., 2002) and NormFinder (Andersen et al., 2004) and were also analyzed with the C q method (Silver et al., 2006).

Gene Expression Variation Analysis
The geNorm algorithm calculates the expression stability value (M) for each gene based on the average pairwise variation (V) of that gene with all the other genes, and determines the number of genes needed for normalization by calculating the pairwise variation (V n /V n+1 ) between two sequential sets containing an increasing number of genes (Vandesompele et al., 2002).
NormFinder is a Microsoft Office Excel ad-in implementing a mathematical model-based variance estimation algorithm to evaluate the overall variation of the genes in the sample as well as variation between the sample subgroups, in order to identify the genes with minimal estimated intragroup as well as intergroup variation (Andersen et al., 2004).
BestKeeper, a Microsoft Office Excel-based tool, ranks the genes based on the average deviation (AD) and the coefficient 1 www.ncbi.nlm.nih.gov of variance (CV) of their C q or C t values across individuals, and estimates the Pearson's correlation coefficient (r) and the probability (P) value for the correlation between each gene and the BestKeeper Index (BKI) represented by the geometric mean of all genes .
The C q method is based on the assumption that the C q value between the two genes remains constant in different individuals when both genes are stably expressed (Silver et al., 2006).
Finally, RefFinder (Xie et al., 2012) is a web-based tool 2 integrating the results of all the four methods for a comprehensive ranking. Unlike the stand-alone versions of the programs, RefFinder only takes quantification cycle values as input, without the possibility to input E values, and efficiency corrected and non-corrected C q and C t values were therefore used in order to assess any differences in the ranking caused by using the efficiency corrected versus non-corrected values.

RNA Quality and Assay Efficiency and Specificity
The average RNA concentration and the optical density ratio of A260/280 were 112.02 ± 39.9 (SD) ng/µl and 1.98 ± 0.06, respectively. The Bioanalyzer electropherogram for each sample showed a clear presence of two peaks expected for the two large rRNA species (18S and 28S), without any evidence of degraded products, and the average RNA Integrity Number was 8.92 ± 0.5. The agarose gel electrophoresis and melting curve analysis revealed that all primer pairs amplified a single PCR product with the expected size, and the sequence analysis and BLAST comparison confirmed the amplicons as the genes targeted by the primers.

Amplification Efficiency and Gene Expression Levels
A summary of the estimated qPCR amplification efficiency and expression level variation for each gene is provided in Supplementary Table S4. The mean amplification efficiency (i.e., the ratio of the number of target gene molecules at the end of a PCR cycle divided by the number of target molecules at the start of the PCR cycle) ranged between 1.66 and 1.92 and the mean C q values between 8.27 and 25.62 (Figure 2). The Rn18s gene showed the most abundant transcript levels (C q of 7.3-9.2), while the Actb gene showed the lowest expression (C q of 23.7-29.6). The largest dispersion of C q values of 5.9 cycle was observed for Actb (Figure 2) and the lowest for Rn18s and Sdha (1.9 cycle for both).

Gene Expression Stability
The M values of gene expression stability among the bank vole samples (n = 21) calculated by geNorm varied slightly depending on whether the second derivative (C q ) or threshold cycle (C t ) data was used and whether 100% efficiency was assumed for 2 www.heartcure.com.au/reffinder Frontiers in Ecology and Evolution | www.frontiersin.org all reactions or not ( Figure 3A). However, all ten genes had an M value below 1.5 (Table 1 and Figure 3A), which is the recommended cut-off for a gene to be considered stably expressed (Vandesompele et al., 2002). Ube2d2a, Ppia, Tbp and Rn18s had the lowest M values based on C q data and Ube2d2a, Ppia, and Rn18s remained among the four most stable genes also based on C t data, regardless of the assumed reaction efficiency ( Figure 3A). The highest M values were consistently obtained for Gapdh and Sdha ( Figure 3A). The V n /V n+1 values were always lower than the recommended cut-off threshold of 0.15 (Vandesompele et al., 2002), indicating that the combination of the two top ranking genes, Ube2d2a and Ppia, is appropriate to use for normalization ( Figure 3B).
Ube2d2a, Ppia, Tbp, and Rn18s had the lowest SVs and therefore ranked as the four most stable genes in the NormFinder analysis (Table 1), including across the sample subgroups defined by gender (males n = 10; females n = 11) and latitude (northern apex n = 7, central n = 8, southern refugial n = 6), respectively (Supplementary Tables S5, S6), whereas Gapdh and Sdha consistently were at the end of the rankings ( Table 1 and Supplementary Tables S5, S6). The grouped analysis by gender suggested Ppia and Ube2d2a or Ppia and Tbp (depending on the RQ estimate used) as the two best genes with minimal intergroup and intragroup expression variation, whereas Tbp and Ube2d2a were consistently the best pair outperforming any single gene when grouped by the latitudinal zone (Supplementary Tables S5, S6).
The AD values of all genes calculated by BestKeeper were below 1 (Figure 3C), indicating low enough variation for the subsequent analysis. Rn18s and Ube2d2a had consistently the lowest AD, followed by Sdha, Tbp, and Ppia in both C q (Supplementary Table S7) and C t analyses ( Supplementary  Table S8), although in a different order, whereas Actb and Ywhaz consistently showed the highest AD ( Figure 3C). Tbp, Sdha and Ube2d2a had the lowest CVs (Supplementary Tables S7,  S8). All genes were significantly correlated (P < 0.01) with the BKI, with the exception of Gapdh and Sdha (Figure 3C and Supplementary Tables S7, S8).
The lowest levels of deviation in the comparative C q analysis were observed for Ube2d2a, Ppia, Tbp, and Rn18s, regardless of whether C q or C t data were used, whereas Gapdh, Sdha and Actb showed consistently the highest deviation ( Table 1 and Supplementary Figure S1; results derived from C t values not shown).
Based on the geometric mean of the rankings given by each method, RefFinder indicated that Ube2d2a and Ppia were the two overall most stable genes and these results did not differ between the analyses using efficiency corrected and noncorrected C q and C t values as input ( Figure 3D). For the noncorrected C q and efficiency corrected C t , the best two genes were followed by Rn18s and Tbp, for the efficiency corrected C q by Tbp and Sdha, and for the non-corrected C t by Actb and Rn18s. In contrast, Tuba1a, Gapdh and Sdha were overall worst ranked ( Figure 3D).

DISCUSSION
It is becoming increasingly evident that RGs commonly used in RT-qPCR experiments may show considerable expression variation (Hernandez-Santana et al., 2016;Jo et al., 2019). Furthermore, when studying gene expression in wild populations there are additional sources of variation due to the natural heterogeneity of the sample (Kozera and Rapacz, 2013;Huang and Agrawal, 2016). In this study we have identified two strong pairs of RGs for the use in RT-qPCR studies of the bank vole, an emerging model species for evolutionary and epidemiological studies .
Although we did not find absolute consensus between the different approaches, gene rankings were similar and three genes were characterized by high expression stability needed to provide the correct normalization, despite the natural heterogeneity of the sample and regardless of the gene expression measure used and whether the efficiency correction was applied or not. Specifically, Ube2d2a, Ppia, and Tbp were identified as the three most stably expressed genes by GeNorm, NormFinder and the C q method, and ranked as three of the top five genes by BestKeeper. A pairwise variation analysis in geNorm further indicated Ube2d2a and Ppia to be a combination sufficient for the normalization. When considering gender and latitudinal zone as sources of variation (Andersen et al., 2004), NormFinder again suggested Ube2d2a, Ppia, Tbp, and Rn18s as the most stable genes, with the best two-gene combinations supported by multiple RQ estimates being those of Ube2d2a and Ppia and Ube2d2a and Tbp, respectively.
Although not traditional housekeeping genes, the function and expression stability of Ube2d2a and Ppia are likely not FIGURE 3 | Gene expression stability analysis. (A) The average expression stability (M) calculated with geNorm, with the more stable genes toward the right-hand end of the graph. For comparison, results based on the relative quantity (RQ) estimated from the second derivative (C q ) data are shown along with those based on the threshold cycle (C t ) data assuming either reaction-specific or 100% qPCR efficiency. (B) Pairwise variation (V n /V n+1 ) calculated with geNorm for two sequential sets containing an increasing number of genes. A value below 0.15 indicates that the addition of another gene has no significant effect in data normalization. (C) BestKeeper ranking by decreasing average deviation (AD) in gene expression values. Pearson correlation coefficient (r) of the expression level of each gene with the BestKeeper Index (the geometric average of all genes) is also shown. For comparison, results based on the C q data are shown along with those based on the C t data. (D) Comprehensive ranking with RefFinder by the geometric mean of the ranking weights assigned to each gene based on the results from geNorm, NormFinder, BestKeeper, and the C q method. Results based on efficiency corrected and non-corrected C q and C t values are shown for comparison. tissue specific to spleen in the bank vole. Both genes show low tissue specificity in other rodents and in humans (Ryffel et al., 1991;Jensen et al., 1995;Wing et al., 1996;Kim et al., 2014), and they play important roles in cell signaling related to post-translational modification and transcription activation, respectively, which are conserved across different cell types and species (Jensen et al., 1995;Wing et al., 1996;Guo et al., 2016).
A recent study of a different mammalian model, yaks, found Ube2d2a and Ppia to be stably expressed across heart, liver, spleen, lung, kidney and skin tissues . From the consistency of our results, it is clear that Ube2d2a and Ppia present higher expression stability in our system than some of the commonly used RGs, i.e., Actb, Gapdh, Sdha, or Tuba1a, and should be considered suitable alternatives. These results The expression stability values were estimated from C q data with BestKeeper and the C q method, and by using efficiency-corrected relative quantity in geNorm and NormFinder.
support that common RGs may show high expression variability in some tissues or experimental settings and other genes may thus present better RGs (Hernandez-Santana et al., 2016;. A possible explanation for a higher expression variation of a housekeeping gene in some tissues is that alternative genes carry the same basal "housekeeping" functions in different tissues (Zhang et al., 2015). The remaining three genes (Rn18s, Rplp0, and Ywhaz) showed good stability values and may also make suitable RGs under certain experimental conditions. For example, because of the high abundancy of ribosomal RNA transcripts in most cells (O'Neil et al., 2013), Rn18s can serve as RG when analyzing levels of highly expressed target genes to ensure similar reaction kinetics (Vandesompele et al., 2002).
Our results are largely concordant with a recent study of Dvořáková et al. (2020), which compared the expression levels of six genes in a small number of bank voles and identified Ppia and Rn18s as low variance RGs. Interestingly, significant differences in the expression levels of the highly expressed bank vole globin genes (mean C q of 12) were only found when data were normalized to Rn18s but not to Ppia (Dvořáková et al., 2020), supporting the importance of similar transcript levels between the RGs and target genes (Vandesompele et al., 2002).
The fact that the genes show stable expression levels in a sample of bank voles from across Europe suggests they may fit as RGs for a broad range of geographical settings. The robustness of the results with respect to the quantification cycle measure (second derivative or threshold) and reaction efficiency assumptions (reaction-specific versus 100%) indicates these genes are also likely to perform well with different protocols and laboratory setups.
The findings of our study are based on data from spleen, which was the available tissue in our broad geographical sampling, but we recognize that it may not be fully representative of every tissue. Despite this limitation, we anticipate that the selected RGs will be applicable to studies addressing a broad range of biological and biomedical problems. The transcription in spleen responds to various physiological and pathological processes, which makes it a suitable tissue for broad range of studies including immunology, toxicology, parasitology and physiology (Guivier et al., 2014;Rossow et al., 2014;Dubois et al., 2018;Dvořáková et al., 2020). The fact that it is a hematopoietic organ in rodents (Brodsky et al., 1966) makes spleen particularly well-suited also for gene expression studies of erythroid and other hematopoietic genes, including globin genes with important ecological and evolutionary functions (Kotlík et al., 2014;Dvořáková et al., 2020).
Taken together, our results highlight the need of a systematic evaluation of appropriate RGs for each study system. Our approach has identified Ube2d2a and Ppia or Ube2d2a and Tbp as the best combinations of RGs in a naturally heterogeneous sample of bank voles while some of the traditional housekeeping genes, such as Actb, Gapdh, Sdha, and Tuba1a, have shown lower stability, and therefore we cannot recommend their use for normalization. It remains unclear to what extent this result applies to other tissues, and researchers planning qPCR geneexpression experiments in the bank vole should make every effort to evaluate the RGs in their particular setting. In this regard, the findings of our study should provide a valuable starting point for such effort.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
LN, SM, and PK designed the study and wrote the manuscript. LN and SM carried out the experiments. LN analyzed the data. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
Petra Šejnohová provided expert technical assistance and Ahmed Gad read the manuscript and offered helpful suggestions.

562065/full#supplementary-material
Supplementary Figure 1 | Gene expression stability analysis using the C q method. Supplementary Table 4 | Summary of variation in the qPCR amplification efficiency (E) and quantification cycle values C q , C t , and efficiency-corrected C t for each gene.

Supplementary
Supplementary Table 5 | Gene ranking based on stability values (SVs) from grouped analysis with NormFinder using efficiency-corrected relative quantities (RQs) derived from the C q data as input.
Supplementary Table 6 | Gene ranking based on stability values (SVs) from grouped analysis with NormFinder using relative quantities (RQs) derived from the C t data and assuming either reaction-specific or 100% efficiencies of qPCR. Table 7 | Expression stability analysis with BestKeeper based on the C q data, including the pairwise correlation of each gene with the BestKeeper Index.

Supplementary
Supplementary Table 8 | Expression stability analysis with BestKeeper based on the C t data, including the pairwise correlation of each gene with the BestKeeper Index.