Selection of reference genes suitable for normalization of qPCR data under abiotic stresses in bioenergy crop Arundo donax L.

Suitable reference gene selection in qRT-PCR is a key pre-requisite to produce reliable data in gene expression analyses. In this study, novel primers for six commonly used reference genes (AC1, TLF, Act2, TUB α, EF-1α and GAPDH) plus two new candidates (pDUF221 and RPN6) were designed and comparatively tested for expression stability under abiotic stresses (osmotic, heavy metal and heat shock) in shoot, root and their combination of Arundo donax L., a raising non-food energy crop. Expression stability rankings from the most to the least stable gene in each condition and in two tissues (young shoots and roots) were generated with geNorm, NormFinder and BestKeeper programs. All programs provided similar rankings and, strikingly, in most cases identified one of the new candidates, RPN6, as the most suitable reference gene. This novel set of reliable references allows to choose either the best combination of reference genes across multiple stress/organ conditions or to select condition-specific genes that can improve the quality of qRT-PCR analysis. This work provides a solid basis for the functional characterization of A. donax, by enabling accurate quantification of the transcriptional responsiveness under a series of common stress conditions of any gene of interest in this promising biomass/bioenergy species.


Results
Gene selection, amplification specificity and efficiency. In this study, we selected eight candidate reference genes from different sources. Four of them, namely TLF (translation factor), Act2 (actin2), Tub α (alpha tubulin) and EF-1 α (elongation factor 1-alpha), have been already used in qRT-PCR experiments in foxtail millet 28 , the species with a sequenced genome which is most closely related to A. donax 23 . For all these genes, the A. donax homologs which were used for primer design shared a nucleotide similarity greater than 91% to foxtail millet references. Published primers for foxtail millet were slightly adapted to giant reed with exception of TLF for which completely new primers were designed (Table 1, Supplementary Table S1). By mining the sorghum transcriptome 29 , we extracted two more candidates: the commonly used GAPDH (glyceraldehyde-3-phosphate dehydrogenase 2) and RPN6 (26 S proteasome non-ATPase regulatory subunit 11). A. donax homologs share respectively 91% and 92% of sequence similarity to their closest sorghum homologs (Supplementary Table S1). The other candidate gene homologous to an Actin gene (AC1) in sorghum (having 89% similarity at nucleotide level) used in the previous study 24 was included in this analysis. Finally, one new gene pDUF221 has also been considered as a candidate based on transcriptomic data analysis in A. donax 24 . To determine the amplification specificity of each primer pair prior to intensive qPCR analysis, electrophoretic analysis on agarose gel and melting curve assessment were carried out for each amplification; a single amplification band with expected size and length for each primer pair was observed on 2% agarose gels, and a single peak was detected through melting curve analysis for all amplifications and no signal for non-template control (Supplementary Figure S1). Moreover, for all primer pairs, amplification efficiencies were calculated based on standard curve assay generated from amplifications with a series of cDNA dilutions and resulted to be between 92.85% and 104.03% with correlation coefficients (R 2 ) ranging from 0.987 to 0.998 (Table 1, Supplementary Figure S2).
Gene expression profile of candidate reference genes. To assess the transcriptional variation and stability of each candidate reference gene, cycle threshold (Cq) value was obtained using qPCR to estimate the expression level of each gene among all the samples (two A. donax tissues for three stress conditions, each with five time points plus one pre-stress control) (Supplementary Table S2). The highest Cq value was detected for AC1 (lowest expression: 28.58 cycles), while the lowest Cq value was measured for GAPDH (highest expression: 16.16 cycles). Mean expression values per gene varied from 25.28 of pDUF221 to 18.84 of GAPDH. To provide a more informative stability index, we calculated also the difference between 75 th and 25 th percentile (ΔP), which was inversely proportional to the spread of the data 30 . Based on this criterion, RPN6 (ΔP = 0.82) was the most stable gene, followed by Act2 (1.07), GAPDH (1.01), TUB α (1.27), pDUF221 (1.31), TLF (1.32), EF-1α (1.48) and AC1 (P = 1.87) (boxes in Fig. 1). This stability ranking among genes is also confirmed by the coefficient SCIEntIFIC REPoRTS | 7: 10719 | DOI:10.1038/s41598-017-11019-0 of variation (CV) that ranged from 2.5% of RPN6 to 6.24% of AC1 (Supplementary Table S3). Comparison of the expression profile of each gene from 0 (control) to 24 h after stress application (Supplementary Figure S3) showed that GAPDH was always the most highly expressed gene in every condition/tissue. The lowest expression levels were mainly associated to pDUF221, with the exception of heat shock stress. Only for this stress, in shoot RPN6 had higher Cq at 3 h and 6 h and, in root, RPN6 and AC1 had higher Cq values for the whole time course (Supplementary Figure S3).

Overall expression stability analyses of candidate reference genes under stress treatments.
All treatments, divided according to tissue, were included in an overall analysis to understand which candidate gene is the most stable across samples. For shoot tissue, the results were heterogeneous, with a prevalence of  F: TCTTGGCTTGCATTCTTGGG  93  100,72 0,998  81,5  24  R: TGGATTGCGAAGGCTGAGTAC   Act2  Actin 2  Unigene057037  F: CGCATACGTGGCACTTGACT  126  92,85  0,987  83,5  28  R: GGGCATCTGAACCTCTCTGC   EF-1α  Elongation factor  1-alpha  Unigene076509  F: TGACTGTGCTGTGCTCATCA  133  97, EF-1α and RPN6 in the top three genes for all the algorithms ( Table 2). In this case, geNorm indicated TLF and TUB α as the most stable genes, while they were ranked as sixth and fourth by NormFinder and as last and fifth by BestKeeper. The least stable genes were instead consistently AC1 and pDUF221 with the only exception of BestKeeper that identified, as mentioned above, TLF as the worst ranking gene. Also, in root EF-1α and RPN6 genes were always in the top three positions (with the exception of Act2 ranked first by BestKeeper), while pDUF221 and AC1 were relegated strictly to the last two ranks by all algorithms. Taking into account both shoot and root organs, RPN6 appeared as the candidate gene with the best scores in the three algorithms followed by EF-1α (first, second and sixth) (Supplementary Table S4). Overall intergroup and intragroup analyses were also carried out for all samples with NormFinder. This analysis graphically showed the variation of each gene with respect to four subgroups: control, osmotic, heavy metal and heat shock stresses (Supplementary Figure S4). Analysis of the best combination of two genes was consistent with single stability value in root, where EF-1α + RPN6 (combination value of 0.105) was the best pair of genes. More surprisingly, in shoot the best pair of genes suggested were RPN6 (second, S = 0.188) and TLF (sixth, S = 0.320) with a combination value of 0.091. It should be noticed that also geNorm algorithm included TLF in the best pair of primers together with TUB α. Considering together shoots and roots, NormFinder suggested the combination between RPN6 and GADPH with a value of 0.115 (the latter classified at the sixth position as a single gene) (Supplementary Figure S4C). These differences with respect to the results of the gene by gene analyses were due to the fact that NormFinder tended to enhance intergroup stability by balancing over/under expression of the two genes to be as close as possible to zero in each subgroup.
Expression stability analyses of candidate reference genes under water/osmotic stress treatment. Water stress was induced by adding 15% of PEG 6000 to the hydroponic solution which decreased the osmotic pressure of the media, consequently inducing osmotic and water limitation stress in the plants 31 without toxic effects in the short term 32 . Reference gene analyses under this simulated drought condition in A. donax shoots showed a high and shared rank for RPN6 (first in NormFinder and BestKeeper, second in geNorm) followed by EF-1α gene (first in geNorm, fourth in NormFinder and third in BestKeeper) ( Table 3). In roots, RPN6 was confirmed as the most stable gene followed by GAPDH that had the best value in geNorm (shared with RPN6) and BestKeeper. Using both organ's samples, RPN6 and GADPH were the best two genes (Supplementary Table S4), resulting in the top three ranks provided by all programs. The lowest stability values resulted from the AC1 and pDUF221 genes with all algorithms and tissues with exception of BestKeeper that ranked pDUF221 at the fourth position in root and TLF and EF-1α as last in the shoot + root analysis.   Expression stability analyses of candidate reference genes under heavy metal stress treatment. A. donax is considered a suitable plant for phytoremediation of contaminated soil. Among the heavy metals that can affect soil quality and reduce plant productivity, cadmium (Cd) is one of the most toxic heavy metals 33 . A. donax seems able, however, to cope with it without physiological adaptation, which is an important feature for phytoremediation 22,34 . When subjected to cadmium toxicity stress, reference gene performances were more heterogeneous respect to osmotic stress: in shoots, RPN6 and TUB α ranked in the top three for all the programs, but in roots TLF was among the best three followed by EF-1α (first, second and fourth) ( Table 4). On the other hand, the least stable gene in shoot was pDUF221, in root was AC1 and in the full set of heavy metal stressed samples were GADPH and AC1 (Table 4). Noteworthy, BestKeeper indicated as most stable candidate pDUF221, which was instead ranked fourth by geNorm and last by NormFinder. Interestingly, analysis of shoot and root together put pDUF221 always in the top three positions, indicating that this gene is the most suitable across different organs for heavy metal stress treatments, but less stable if organs are taken separately (Supplementary Table 4).
Other general candidate genes to be used in studies encompassing both root and shoot could alternatively be the common Act2 (ranked fourth, second and first) and TLF (ranked third, first and third).
Expression stability analyses of candidate reference genes under heat shock stress treatment. Heat shock is an important factor that affects plant physiology and growth 35 . The capability of A. donax to survive in warm environments is an interesting trait of this plant 16 , which can be relevant to forecast its productivity as heat spells become more frequent. Gene stability in heat condition (42 °C) identified the RPN6 gene as the most stable reference gene in shoot (third, first and first for geNorm, NormFinder and BestKeeper, respectively), root (first, first and second) and also considering both organs together (always first) ( Table 5). Other suitable reference genes for this stress condition were EF-1α (first, second, second) and TUB α (second, third, fourth) in shoot, EF-1α (first, second, first) and Act2 (third, third and fourth) in root and TLF (first, second, fourth) in overall heat shock (Supplementary Table S4). Again, the genes with overall lower ranking across algorithms were pDUF221 and AC1 in all the conditions.
Best reference gene number identification. geNorm Excel add-in is a useful tool for calculation of the best number of genes that should be used together in a relative qRT-PCR experiment. Therefore, it has been used in this study to predict the optimal number of reference genes to be used in each stress experiment (osmotic, heavy metal and heat shock). As expected, the number of references depended on the experimental settings:   considering shoot and root together geNorm indicated three genes as the most suitable for single stresses and four genes considering all stresses ( Fig. 2A). If root and shoot were considered separately, all the values for single stress dropped below the suggested threshold of 0.15 (Fig. 2B,C). With one exception in shoot, by grouping the three stresses, the value was just above the threshold limit when two candidate genes were used while it laid on the threshold with three.

Validation of reference genes by expression pattern analyses of target genes under stress treatments.
To demonstrate the reliability of the newly analysed reference genes in A. donax for stress treatments, DREB2A (Dehydration-Responsive Element Binding Protein 2 A) which is one of the key genes triggering the response to both drought and heat shock 26,36,37 and isoprene synthase gene (IspS), which is differentially expressed in heavy metal (Cd) stress treatment in Arundo donax 27 , were chosen to validate these reference genes for qPCR internal normalization. We compared the expression profiles of DREB2A in shoot under osmotic stress and heat shock treatments using the two most stable genes (RPN6 and EF-1 α), their combination (RPN6 + EF-1 α) and the least stable genes (AC1 and pDUF221) as reference genes based on stability ranking (Tables 3 and 5, respectively). Fold change of DREB2A was calculated with the comparative Ct method 3 . The expression pattern was consistent with a two-fold increase of DREB2A expression at 6 h and 11 h and about four-fold increase at 24 h normalized with RPN6, EF-1α and their combination (Fig. 3A). On the other hand, the pattern obtained using AC1 as reference displayed an upwards-shifted trend, with a twofold increase since the beginning and a six-fold increase from 6 h on. When pDUF221 was used as a reference gene, the expression pattern decreased at 3 h, 6 h, and 11 h to finally grow again at 24 h. Under heat shock treatment, the expression levels of DREB2A normalized with RPN6, EF-1α and their combination were extremely consistent at all time points: after a 12-fold increase at 1 h30′, DREB2A expression gradually decreased to four-fold at 3 h and two-fold at 6 h. By contrast, DREB2A expression levels were found to be two-fold higher at 1 h30′, 3 h and 6 h using AC1 as reference gene and at least two-fold lower at all time points using pDUF221 compared to when using RPN6, EF-1α and their combination (Fig. 3B). Furthermore, the expression levels of IspS under cadmium treatment were highly consistent at all time points with a 1.8-fold increase at 1 h30′ and a four-fold increase at 3 h using RPN6, TUB α and their combination as reference genes, respectively. Meanwhile, IspS expression levels were shown to decrease by 30% at all time points when normalized with TLF compared with RPN6, TUB α and their combination. When the expression levels of IspS were normalization with pDUF221, we observed a dramatic increase. A 2.8-fold increased was observed at 1 h30′, followed by a continuous increase from 6 h till 24 h, instead of the decreased expression at 11 h and the further increase at 24 h observed when normalizing with RPN6, TUB α and their combination (Fig. 3C).

Discussion
Real Time PCR is a powerful technique for gene expression profiling and functional characterization of genes. Its accuracy, however, is critically dependent from the choice of reference genes whose expression is strictly proportional to the total mRNA amount in the samples to be quantified 38 . Therefore, ideal reference genes should have the same expression levels (measured in terms of Cq value) among different conditions, tissues, developmental stages and crop varieties 7 . In reality, however, gene stability depends greatly on the plant-environment system so that different stress conditions, organs or cultivars can bring undesirable variations in expression of reference genes selected without specific validation, possibly leading also to erroneous quantification 8 . For this reason, a full set of possible reference genes should be developed and systematically tested for each species whenever the experimental conditions change. In this study, eight candidate reference genes were selected for A. donax L., an emerging non-food energy crop 39 : AC1, pDUF221, TLF, Act2, TUB α, EF-1α, RPN6 and GAPDH. The aim of this study was not only to evaluate the overall performance of these genes as references for qRT-PCR data analysis, but also to provide a detailed indication on which gene is more suitable for specific organ/stress combinations. Identifying stable reference genes is not a trivial task, as it encompasses the choice of the candidates based on literature searches, the identification of the homologs from a non-model species, the design of the primers and their thorough verification in the species of interest, in this case A. donax. Any of these steps can fail, resulting in a considerable waste of time, effort and resources. As long as the conditions used in future studies will be the same as those used here, therefore, the relevance of our results is that they provide for A. donax reliable references without the need for time-consuming testing/optimiziation. Given its biological features and economic relevance 24,25,40 , the three kinds of abiotic stress applied in this study (osmotic, heavy metal and heat shock) are among the most interesting ones to characterize the functional bases of A. donax tolerance to adverse environmental conditions. Three algorithms were applied to the Cq values to measure the stability of the candidate genes: geNorm, NormFinder and BestKeeper. Our results show overall congruence on the stability ranking of the reference genes generated from different algorithms, which indicates the good performance and reliability of the methods. We also divided the analysis into stress and organ subsets to identify for each subset of samples the best performing reference genes. In general, gene RPN6 was always ranked among the top three genes in both organs under osmotic stress, heat shock and overall analysis, which makes it a highly suitable reference (Table 2). Nevertheless, RPN6 showed under heavy metal stress, in root and with the complete dataset slightly higher variability, especially when the NormFinder and BestKeeper algorithms were used. On the other hand, geNorm classified RPN6 as the best reference gene also in root, making the choice difficult in this context. These results are very interesting, giving the fact that until now, this gene was considered as possible reference gene in quantitative RT-PCR only in Arabidopsis 41-43 but, to our knowledge, never in monocot species. Its putative orthologue in Arabidopsis (AT1G29150) encodes a protein constituting a lid subunit of the 26 S proteasome, which is involved in the ubiquitin-proteasome system (UPS) for degradation of misfolded proteins and for stress response 44 . Other subunits of this large protease complex have been previously suggested as potential source of new and more stable reference genes in Arabidopsis 45 . Worth of notice, in Arabidopsis RPN6 is classified as a cadmium responsive gene 46 , which could explain the variation we observed in response to this stress. We notice, however, that RPN6 still remained more stable than some other commonly used reference gene such as EF-1α or AC1 (Table 4, Supplementary Table S4), indicating a pretty limited transcriptional resposiveness of RPN6 to cadmium in A. donax as compared to A. thaliana. This fact is possibly due to the high resistance of A. donax to cadmium treatment as demonstrated by Papazoglou and colleagues 34 . A useful reference gene across different stresses can be the commonly used EF-1α, which ranked always in the top three positions. From our analysis, this gene often scored better when the organs were considered separately, especially in osmotic and heat shock stresses (Tables 3 and 5), suggesting moderate tissue-specific variability. The EF-1α gene from Setaria italica (L.) P. Beauvois (foxtail millet) was suggested as the best internal control gene under drought and salt condition 28 . Here, we confirm the good stability of gene EF-1α also in A. donax, a close relative of S. italica, but mainly when used in single tissues. Considering both tissues together, the Act2 gene performs better in each stress condition, especially under heavy metal treatment (Supplementary Table S4). GAPDH is, instead, suitable in osmotic stress, where it resulted the best among commonly used reference genes (Table 3). This result is consistent with previous studies showing that GAPDH is one of the best references in plants under drought and salt stress, but not in other conditions [47][48][49] . Finally, we suggest to avoid using two of the candidate genes tested due to their instability: pDUF221 and AC1. In the case of pDUF221 this result is likely due to fact that the transcriptomics analysis used for its selection encompassed a single time point only from osmotic stress, while here we tested a 24 h time series and multiple stresses. As for the number of reference genes to be used, geNorm analyses showed some interesting trends based on the different dataset employed. In general, the higher the number of organs/conditions, the higher tended to be the number of genes required for a reliable quantification. In particular, the full dataset (all the treatments + both organs) dropped below the suggested 0.15 threshold only with four reference genes, while in roots it estimated that only two genes were needed ( Fig. 2A,C). This is significant because it implies a differential regulation between organs that becomes difficult to quantify in case of studies aimed at the comparison of target genes of aerial and underground tissues. Another interesting point is the marginal decrease of the pairwise variation when considering the combination of all stresses in shoot with two or three reference genes (Fig. 2B). Given the relatively high cost/benefit ratio in the use of an additional reference gene, and considering that 0.15 is an arbitrary threshold 10 , we suggest that use of two internal controls can be a good compromise for both root and shoot when interested in the organ-specific comparison of target expression levels under osmotic, heavy metal and heat stress in A. donax. DREB2A is a well-known transcription factor, which is associated with drought, salt and heat responses. The DREB2A protein interacts with a cis-acting dehydration-responsive element and activates a downstream cascade of drought and heat-responsive genes, thus providing a better tolerance for plants to these stresses 26,36,37 . The consistent expression pattern obtained for DREB2A under osmotic and heat shock treatment normalized with the two best reference genes and highly variable pattern with the two worst ones further confirmed the reliability of the ranking of these reference genes. For the cadmium stress we used the IspS gene, which we recently characterized in A. donax and found to be transcriptionally upregulated by cadmium treatment 27 . The results of the validation indicate that, in the case of heavy metal, the performance of the reference genes for normalization between the best and least stable genes does not differ as much as in the cases of those for PEG and heat shock treatments. Given the relatively low cadmium responsiveness of IspS, however, it is possible that the differences in normalization among reference genes may result bigger for genes with stronger upregulation in response to cadmium.
In summary, this study provides a wide view of the reference genes that can be used or should be avoided in A. donax under specific abiotic stresses and in specific organs, making an important step forward towards the reliable and accurate gene expression quantification in this species. In addition, this study emphasizes further that normalization with reference genes rigorously validated before use for any new experimental design is essential. Moreover, thanks to analysis of related species transcriptomes, a new stable gene (RPN6) has been successfully used for relative quantification, showing that a deeper comparative analysis of plant transcriptomes can unveil additional candidates for a more precise and reliable qRT-PCR analysis. We expect that RPN6 could find application in additional monocotyledonous species.

Methods
Plant materials and stress treatments. Cohorts of A. donax cuttings (Sesto Fiorentino, Florence, Italy 43°49′01.8′N 11°11′57.0′E) were used in this study. The plant growing condition and procedure for stress treatments were the same as those described previously 24 . For stress treatments, plants at the 5-leaf stage were transferred from hydroponic solution to fresh one supplemented with 15% PEG 6000 (osmotic stress), 500 µM CdSO 4 (heavy metal stress) or pre-warmed at 42 °C (heat stress). The entire shoots and roots (treated and untreated) were independently collected at all different time points (0 h, 1 h30', 3 h, 6 h, 11 h and 24 h), immediately frozen in liquid nitrogen, and then stored at −80 °C till use. Three biological replicates were applied for all the treatments at every sampling time point.
SCIEntIFIC REPoRTS | 7: 10719 | DOI:10.1038/s41598-017-11019-0 Candidate reference gene selection, PCR primer design. Among the sequenced genomes deposited in Phytozome, Sorghum bicolor and foxtail millet (Setaria italica L.) are the two species phylogenetically most closely related to A. donax 23 . A.donax homologs of four common housekeeping genes from foxtail millet TLF (Terminal Flower-like), Act 2 (Actin 2), Tub α (Tubulin alpha), and EF-1 α (Elongation Factor 1 alpha) 28 were selected by BLASTN searches against the giant reed reference and water-stress transcriptomes 23,24 . Analogously, two additional genes were identified in the A. donax transcriptomes based on stable expression across tissues (RPN6 29 ) or because already used as a reference gene in sorghum (GAPDH). Finally, the last two genes were selected directly from the giant reed transcriptome (AC1 and pDUF221 24 ) based on their low coefficient of variation (CV) across organs/water stress conditions. The AC1 candidate had been previously chosen (Fu et al. 2016), because among the different A. donax transcripts coding for actin, it was the one with the lowest CV in the PEG-treated transcriptomes. The pDUF221 candidate was chosen exclusively based on the CV from the PEG-treated transcriptomes, where it ranked 5 th among the most stable transcripts for this stress (Supplementary Figure S5).
Total RNA isolation and cDNA synthesis. Total RNA was isolated with the Spectrum Plant Total RNA Extraction Kit (Sigma) for shoots and the Rneasy ® Plant Mini Kit (Qiagen) for roots, respectively. To assure complete absence of genomic DNA contaminations, extracted total RNA was treated with DNase I (Sigma-Aldrich) and checked on 1% agarose gel for integrity control. Concentration and quality of each sample were measured spectrophotometrically through the OD 260 /OD 280 absorption ratio. First strand cDNA was reversed transcribed from 1 µg of total RNA primed with oligo-dT in a total reaction mixture of 20 µL using SuperScript ® III Reverse Transcriptase (Life Technologies) according to the manufacture's instruction.

PCR amplification specificity and Quantitative Real Time PCR analyses.
To assess the amplification specificity of each primer pair prior to qPCR analysis, PCR amplification was performed in a total volume of 10 µL containing 1 µl of 6-fold diluted cDNA (8 ng Figure S1B). A standard curve of qPCR reaction was generated from five points (four points for pDUF221gene and target gene DREB2A) of a 6-fold dilution series (10-fold dilution for pDUF221 and TLF). The slope (S) of the standard curve was used to calculate the amplification efficiency (E) of each primer pair as follows: E = 10 (−1/S) (Supplementary Figure S2). Three technical replicates were used for each sample and every plate contained one No Template Control (NTC) well for each primer pair used. In order to compare different plates, in Bio-Rad CFX Manager software the baseline threshold was set at 329.82 and one control sample was used in every plate to check for Cq congruency. For the validation of the reference genes under cadmium stress, as IspS is under circadian regulation 50 , at each time point the relative expression was normalized with the untreated control at the same time.
Data analyses with geNorm, NormFinder and BestKeeper. The analyses were conducted on six different datasets that comprise each a single combination of organ (shoots or roots) and indpendent stress treatment (osmotic, heavy metal or heat shock); in addition, one analysis was performed on the combined dataset comprising both organs in all stress conditions. This is necessary because the reference gene may vary depending on the experimental settings 7 . Three different Excel-based algorithms have been applied for data analysis. geNorm v3.5 10 and NormFinder v0.953 11 require relative input data, so the Cq values were converted with the formula 2 −ΔCt where ΔCt is the difference of each Cq value minus the lowest Cq value (highest expression level). BestKeeper 12 instead, uses raw Cq values. geNorm calculates stability value (M) based on the average pairwise comparison with a stepwise exclusion of the highest M value (least stable gene). Further, geNorm calculates the number of genes needed for a reliable normalization considering the pairwise variation (V n /V n+1 ) between sequential normalization factors, NF n and NF n+1 . This number is optimal when the addition of one more reference gene does not significantly contribute to the variation of the normalization factor (NF n+1 ) or, as suggested, the value drops below 0.15. NormFinder uses an ANOVA-based algorithm to estimate intra-and inter-group variation for a given set of experiments, providing a rank where the most stable gene is the one with lowest stability (S) value. Moreover, NormFinder provides the best gene pair combination that minimizes the expression differences among subgroups, if subgroups are set. BestKeeper, differently from the other algorithms, does not provide a direct rank list but calculates standard deviation (SD [ ± Cq]) and coefficient of variation (CV [%Cq]) for each gene. We sorted the CV values to rank the genes from most stable (lowest CV value) to least stable (highest CV value).

Data Availability
All data analysed during this study are included in this article (and its Supplementary Information files).