Selection and Validation of Reference Genes for Accurate RT-qPCR Data Normalization in Coffea spp. under a Climate Changes Context of Interacting Elevated [CO2] and Temperature

World coffee production has faced increasing challenges associated with ongoing climatic changes. Several studies, which have been almost exclusively based on temperature increase, have predicted extensive reductions (higher than half by 2,050) of actual coffee cropped areas. However, recent studies showed that elevated [CO2] can strongly mitigate the negative impacts of heat stress at the physiological and biochemical levels in coffee leaves. In addition, it has also been shown that coffee genotypes can successfully cope with temperatures above what has been traditionally accepted. Altogether, this information suggests that the real impact of climate changes on coffee growth and production could be significantly lower than previously estimated. Gene expression studies are an important tool to unravel crop acclimation ability, demanding the use of adequate reference genes. We have examined the transcript stability of 10 candidate reference genes to normalize RT-qPCR expression studies using a set of 24 cDNAs from leaves of three coffee genotypes (CL153, Icatu, and IPR108), grown under 380 or 700 μL CO2 L−1, and submitted to increasing temperatures from 25/20°C (day/night) to 42/34°C. Samples were analyzed according to genotype, [CO2], temperature, multiple stress interaction ([CO2], temperature) and total stress interaction (genotype, [CO2], and temperature). The transcript stability of each gene was assessed through a multiple analytical approach combining the Coeficient of Variation method and three algorithms (geNorm, BestKeeper, NormFinder). The transcript stability varied according to the type of stress for most genes, but the consensus ranking obtained with RefFinder, classified MDH as the gene with the highest mRNA stability to a global use, followed by ACT and S15, whereas α-TUB and CYCL showed the least stable mRNA contents. Using the coffee expression profiles of the gene encoding the large-subunit of ribulose-1,5-bisphosphate carboxylase/oxygenase (RLS), results from the in silico aggregation and experimental validation of the best number of reference genes showed that two reference genes are adequate to normalize RT-qPCR data. Altogether, this work highlights the importance of an adequate selection of reference genes for each single or combined experimental condition and constitutes the basis to accurately study molecular responses of Coffea spp. in a context of climate changes and global warming.

World coffee production has faced increasing challenges associated with ongoing climatic changes. Several studies, which have been almost exclusively based on temperature increase, have predicted extensive reductions (higher than half by 2,050) of actual coffee cropped areas. However, recent studies showed that elevated [CO 2 ] can strongly mitigate the negative impacts of heat stress at the physiological and biochemical levels in coffee leaves. In addition, it has also been shown that coffee genotypes can successfully cope with temperatures above what has been traditionally accepted. Altogether, this information suggests that the real impact of climate changes on coffee growth and production could be significantly lower than previously estimated. Gene expression studies are an important tool to unravel crop acclimation ability, demanding the use of adequate reference genes. We have examined the transcript stability of 10 candidate reference genes to normalize RT-qPCR expression studies using a set of 24 cDNAs from leaves of three coffee genotypes (CL153, Icatu, and IPR108), grown under 380 or 700 µL CO 2 L −1 , and submitted to increasing temperatures from 25/20 • C (day/night) to 42/34 • C. Samples were analyzed according to genotype, [CO 2 ], temperature, multiple stress interaction ([CO 2 ], temperature) and total stress interaction (genotype, [CO 2 ], and temperature). The transcript stability of each gene was assessed through a multiple analytical approach combining the Coeficient of Variation method and three algorithms (geNorm, BestKeeper, NormFinder). The transcript stability varied according to the type of stress for most genes, but the consensus ranking obtained with RefFinder, classified MDH as the gene with the highest mRNA stability to a global INTRODUCTION Impacts from recent climate-related extremes, such as heat and cold waves, droughts and strong rainfall events, reveal remarkable vulnerability of agricultural systems. Estimated future climate changes are expected to further amplify the existing climate-related risks and create new ones (IPCC, 2013). Within this context, it is important to assess crop ability to adapt their vital processes at a speed compatible with these environmental changes, allowing the selection and breeding of elite genotypes.
Coffee is one of the most important agricultural traded commodities, growing in more than 60 tropical countries. It is estimated that ca. 25 million farmers produce coffee in over 1 million ha, with a majority of smallholders whose livelihoods depend on this crop (Waller et al., 2007). In the last years, worldwide crop yields were above 8 million tons of green coffee beans (ICO, 2014a), generating an income of ca. US$ 173,000 million (ICO, 2014b), and involving approximately 100 million people considering the entire chain of value of coffee (Bunn et al., 2015).
Coffee production and quality are highly dependent on a regular sequence of climate events; temperature and water availability are considered the most important limiting environmental variables to this crop (DaMatta and Ramalho, 2006). Also, several modeling studies on the global impact of climate changes, mostly focused on increased air temperatures, have predicted reductions of suitable areas for Arabica coffee (up to ca. 50%) by 2,050 (Bunn et al., 2015;Magrach and Ghazoul, 2015), compromising the livelihoods of millions of small householders.
Under field conditions, the superimposition of environmental events is the most common situation, as is expected to be the case of elevated [CO 2 ] and enhanced temperatures predicted along this century. In fact, depending on anthropogenic greenhouse gas emission scenarios, [CO 2 ] might reach between 421 and 936 µL CO 2 L −1 by 2,100, accompanied by a global warming up to between 2.6 and 4.8 • C relative to 1986-2005(IPCC, 2013. However, it was recently shown that increased [CO 2 ] in the atmosphere strengthens the coffee plant by reinforcing its Abbreviations: RT-qPCR, quantitative real-time polymerase chain reaction. photosynthetic performance (Ramalho et al., 2013;Ghini et al., 2015;DaMatta et al., 2016). Furthermore, enhanced [CO 2 ] has a crucial role in the mitigation of heat impact, and, therefore, in the resilience of coffee plant to supra-optimal temperatures, with positive repercussions ranging from mineral nutrition (Martins et al., 2014) to the triggering of defense mechanisms and altered gene expression (Martins et al., 2014;Rodrigues et al., 2016). Therefore, the catastrophic predictions on the future of the coffee crop, that are based almost exclusively on temperature drift , should be reconsidered. In addition, results from Rodrigues et al. (2016) strongly suggest the need to study the single and superimposed effects of elevated [CO 2 ] and supra-optimal temperatures at all plant levels, from morphology up to the molecular assessment.
The release of Coffea sp. expressed sequence tag (EST) databases (e.g., Poncet et al., 2006;Vieira et al., 2006;Mondego et al., 2011) and Coffea canephora genome (http://coffeegenome.org/) has greatly prompted the study of genes involved in important agronomic and stress tolerance traits, making marker-assisted selection a straight forward approach. Gene expression analysis is an important tool to elucidate the complex regulatory networks of the genetic, signaling, and metabolic mechanisms that underlie plant-environmental interactions (Mallona et al., 2010). Monitoring differential gene expression and validating high-throughput RNA sequencing (RNA-seq) data is ideally achieved through quantitative real-time PCR (RT-qPCR) analysis. However, accuracy and reliability of RT-qPCR relies on the normalization of gene expression data (Artico et al., 2010;Die et al., 2010). To avoid severe pitfalls in data analysis and interpretation, this implies the selection and systematic validation of suitable reference genes to be used as internal controls. In fact, expression stability should be validated for each particular plant tissue/organ, cell, and experimental conditions, particularly when involving environmental stressful conditions (Die et al., 2010;El Kelish et al., 2014;Imai et al., 2014;da Costa et al., 2015). In addition, the stability of reference gene transcripts is also species-dependent (Andersen et al., 2004;Gutierrez et al., 2008a,b), as in the case of coffee (Cruz et al., 2009;Goulao et al., 2012). A number of housekeeping genes (e.g., b-actin, elongation factor1a, 18S ribosomal RNA, and polyubiquitin) involved in general cell metabolism pathways are widely used to calibrate RT-qPCR studies in biological systems (Willems et al., 2008;Goulao et al., 2012;Imai et al., 2014;da Costa et al., 2015;Llanos et al., 2015). However, their expression levels may vary with the experimental conditions (Petriccione et al., 2015), and their use as internal reference genes should be taken with caution (Gutierrez et al., 2008a,b). For instance, GAPDH has been indicated as one of the most stable genes during single abiotic stresses (Tian et al., 2015), whereas it appeared to be the least stable gene during the plant development . In addition, the speciesdependent stability of reference gene transcripts was observed, for example, in TUB-A, which is the most stable reference gene during celery development (Li et al., 2016) in contrast with what has been observed during the development of cherry (Ye et al., 2015).
In a context of advancing our knowledge regarding the plant responses to estimated climate changes and global warming, this work examined, for the first time, the transcript stability of candidate reference genes to accurately perform the normalization of RT-qPCR data from expression studies of Coffea spp. plants exposed to conditions mimicking predicted future environmental conditions. For that, single and combined impacts of elevated [CO 2 ] and heat were considered.

Plant Material and Growth Conditions
Plant materials and experimental design was previously described in detail . Briefly, three widely cropped coffee genotypes from the two main producing species were used: Coffea arabica L. cv. Icatu (an introgressed variety from C. canephora Pierre ex A. Froehner), C. arabica L. cv. IPR108, and C. canephora cv. Conilon Clone 153 (CL153). Potted plants of 1.5 years in age were transferred from a greenhouse (ambient [CO 2 ]) into walk-in growth chambers (EHHF 10,000, ARALAB, Portugal) differing in air [CO 2 ] supply: 380 µL CO 2 L −1 (380-plants) or 700 µL CO 2 L −1 (700-plants). Both groups of plants were then grown for 10 months in 28 L pots for ca. 10 months under controlled environmental conditions of temperature (25/20 • C, day/night), RH (75%), irradiance (ca. 700-800 µmol m −2 s −1 ), photoperiod (12 h), without limitations of nutrients (Ramalho et al., 2013), water and space for root growth, the latter evaluated by visual examination at the end of the experiment, after removing the plants from their pots.
After that 10 months period, air temperature was gradually increased at a rate of 0.5 • C day −1 , from 25/20 • C up to 42/34 • C, with a 7 day stabilization at 31/25, 37/30 and 42/34 • C to perform several evaluations (Martins et al., 2014Rodrigues et al., 2016). Leaf material was collected from newly matured leaves, of both plagiotropic and orthotropic branches of the upper (illuminated) part of the plant, from five to eight plants per genotype and immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction.

Total RNA Isolation and cDNA Synthesis
Total RNA was isolated and quantified as described in Goulao et al. (2012) for each of the following 24 conditions: three genotypes (CL153, Icatu, IPR108), four temperature regimes (25/20 • C, 31/25 • C, 37/30 • C, 42/34 • C), and two [CO 2 ] (380, 700 µL CO 2 L −1 ), using the RNeasy Plant Mini Kit (Qiagen, Germany) according to the manufacturer's protocol. In order to avoid eventual genomic DNA contaminations, RNA samples were digested with Ambion R DNA-free TM DNase (Ambion, USA). RNA integrity was verified in 1.5% agarose-TAE gel electrophoresis containing GelRed Nucleic Acid Gel Stain (Biotium, USA). All RNA samples were individually analyzed for the possible presence of DNA contamination by standard PCR reactions (35 cycles) using primers designed for ubiquitin (UBQ) gene (Table 1), in the absence of cDNA synthesis. Total RNA concentration and purity were further verified through BioDrop Cuvette (BioDrop, UK) measurements to guarantee the same amount of starting material in subsequent cDNA synthesis. RNA concentrations ranged between 333 and 1,056 ng/µl, with OD 260:280 ratios always above 1.98.
One microgram of DNA-free total RNA was used to synthesize first-strand cDNAs using oligo-(dT) 18 primer and the SuperScript II first-strand synthesis system (Invitrogen, USA).

Quantitative Real-Time Polymerase Chain Reaction Conditions
RT-qPCR reactions were performed in 96-well plates using iQ TM SYBR R Green Supermix (BioRad, USA). Reactions were prepared in volumes of 25 µl containing 150 ng of cDNA and 3 µM each primer, in 1× iQ TM SYBR R Green Supermix. Reactions were run on the iQ TM 5 Real-Time Detection System (BioRad, USA) using the following parameters: hot start activation of the TaqDNA polymerase at 95 • C for 10 min, followed by 40 cycles of denaturation 95 • C for 15 s; annealing at 60 • C for 30 s; elongation at 72 • C for 30 s; and plate read. To verify the specificity of each amplification, and the absence of primer dimers, dissociation curves were obtained for each amplicon at the end of the PCR run, by continuous fluorescence measurement from 55 • to 95 • C, with sequential steps of 0.5 • C for 15 s. Three technical replicates were used for each biological replicate and the mean Ct was used for data analyses. The full sample set was always included in each technical replicate to exclude any artifacts consequential of between-run variations. No signals were detected in non-template controls run in parallel for each primer set. Two negative controls were included for each primer pair, in which cDNA was replaced by water or total RNA.

Calculation of PCR Efficiencies
Five-fold serial dilutions (1:1-1:625) of pooled cDNAs, which included equal-molar quantities of all samples independently for each genotype, were quantified in triplicates to generate standard curves for each primer pair. Based on the slopes of the standard curves, the PCR efficiency of each gene, for each genotype, was determined from the respective logarithm of the cDNA dilution, plotted against the mean threshold cycle (Ct) values. The reaction efficiency was calculated using the equation: Costa et al., 2015), where E is the efficiency, in percentage, and slope is the gradient of the best-fit line in the linear regression.

Selection of Reference Genes
The global transcript variability of each candidate reference gene was analyzed through statistical parameters using Box Plot. The variation of gene expression was calculated by means of coefficient of variation (CV), where CV = (standard deviation mean) × 100. These analyzes were performed using the statistical program R Core Team (2016). In order to select the best reference genes for the experimental conditions, expression stability was analyzed using geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004) and BestKeeper (Pfaffl et al., 2004). The final rank of reference genes was determined with the RefFinder program (Xie et al., 2012), a webbased user-friendly comprehensive tool that integrates geNorm, Normfinder, BestKeeper, and the comparative Ct method.

Determination of the Minimum Number of Reference Genes for Data Normalization
The pairwise variations (V-values) were calculated using the final ranking determined by RefFinder, through pairwise variation (V n /V n + 1 ) of two consecutively serial log 2 -transformed normalization factor (NF) ratios, based on stepwise addition of the subsequent more stable reference gene (NF n relative to NF n+1 reference genes; Vandesompele et al., 2002). The SD of the array generated by the log 2 -transformed NF ratios was calculated for each NF combination (V n /V n+1 ). V n /V n+1 were then plotted, representing changes in expression stability of NF. The optimal number of reference genes to exact standards was determined by calculating values as a pair of variation (Vn/Vn + 1) between two normalization factors consecutively classified (NC) after the gradual addition of more subsequent stable reference gene (NFN and NFN + 1), as proposed by Vandesompele et al. (2002), included in GeNorm package. A cut off value <0.15 was adopted (da Costa et al., 2015;Niu et al., 2015;Petriccione et al., 2015;Yang et al., 2015).

Validation of Reference Gene Analysis
The evaluation of the transcript levels from the RLS gene, encoding the ribulose 1,5-bisphosphate carboxylase/oxygenase (RuBisCO) large subunit, along the imposed temperatures (25/20 • C, 31/25 • C, 37/30 • C, 42/34 • C), and under both growth [CO 2 ] (380, 700 µL CO 2 L −1 ) was analyzed by RT-qPCR, using the most stable vs. the least stable genes (using one, two or three). The RT-qPCR amplification conditions were the same as described above. The relative expression data and statistical analysis were performed using the REST 2,009 software developed by Pfaffl (Technical University Munich) and Qiagen (Germany). Rest analysis takes into account PCR efficiencies and normalization gene, using the test hypothesis P (H1), i.e., the probability of difference between target and control groups be by chance. P (H1) performs at least 2,000 times of random reallocations of target samples and controls between the groups. Statistical differences were considered significant when P < 0.001 ( * * * ), P < 0.01 ( * * ), and P < 0.05 ( * ).

RESULTS AND DISCUSSION
A total of 10 candidate genes (CG) to obtain an accurate validation of RT-qPCR data were assessed considering the individual comparison groups, i.e., genotype, [CO 2 ] and temperature, multiple stress interaction ([CO 2 ] and temperature) and total stress interaction (genotype, [CO 2 ] and temperature).
To determine the specificity of the primer pairs used in this study, melting curve analysis was performed after 40 cycles of amplification. The single peak obtained confirmed the specificity of the amplicon. No signal was detected in the negative controls of all tested primers. The standard curve method using a pool of all cDNA samples was performed to calculate RT-qPCR efficiency (E) and the correlation coefficient (R 2 ) of each primer pair ( Table 2). The raw quantification cycle (Ct) values ranged from 21.2 to 38.2 across all of the tested samples, with their minimal and maximal means falling between 24.1 ± 1.8 and 34.1 ± 2.1 (Figure 1). The results showed that GAPDH and α-TUB were the most and least abundantly expressed genes with mean Ct of 24.4 and 30.7, respectively. EF-1A showed the most variable levels of expression. Collectively, these results indicate that the transcript levels of the several candidate genes (CG) varied across different experimental samples confirming the need to screen the most appropriate reference genes for each experimental condition to accurately normalize RT-qPCR gene expression data. This result is common in plants given that most reference genes are involved in functions other than basal cell metabolism (Goulao et al., 2012;da Costa et al., 2015). According to the Coefficient of Variation (CV) method, the most stable genes were ACT (for genotype and total stress, CV = 5.91 in both cases), UBQ2 (for temperature, CV = 5.50 and [CO 2 ], CV = 5.71), and S15 (for multiple stress, CV = 5.58; Table 3). EF-1A was ranked as the least stable gene for genotype (CV = 8.49), temperature (CV = 8.34) and overall stress (CV = 8.49), ranking in the 6th place for [CO 2 ] and multiple stress. DNAJ and TUB were the less stable genes for [CO 2 ] (CV = 9.68) and multiple stress treatments (7.95). Because the CV method alone is not sufficiently accurate to estimate the stability of the mRNA levels (Goulao et al., 2012), a complementary analysis was performed based on three algorithms, geNorm, NormFinder, and BestKeeper. According to the geNorm algorithm, the average expression stability (M-value) of all CG was below the default limit of 1.5 for genotype, temperature, [CO 2 ] and their respective interactions, therefore indicating a considerable high stability level (Table 4). GAPDH was consistently classified as the most stable gene for the five comparison groups while CYCL was the least stable gene for all groups with the exception of genotype, where MDH was the least sable gene. On the other hand, the ranking obtained with NormFinder was clearly dependent on the comparison group: MDH was the most stable for genotype, UBQ2 for temperature, GAPDH for [CO 2 ] and ELFA for multiple and total stress groups. Such differences are due probably to different stability measures, i.e., while geNorm calculates gene expression stability based on the average pairwise expression ratio, NormFinder estimates the overall expression variation to provide individual stability values (Kozera and Rapacz, 2013). Nevertheless, similarly to geNorm, NormFinder also ranked CYCL as the least stable gene for all groups, with the exception for genotype (where α-TUB was the least stable CG). In order to consolidate the data obtained with geNorm and NormFinder, we have included in the analysis the BestKeeper algorithm, which calculates the expression stabilities of the CG in accordance with SD and CV values. The calculated results classified ACT as the most stable gene for the genotype (together with MDH) and S15 for the [CO 2 ] group. For temperature, multiple stress and overall stress groups the most stable gene was MDH. α-TUB, EF-1A, and CYCL were ranked as the least stable genes, always below the 6th place and usually in the three last positions. Despite the differences between the ranking obtained with the CV method and the three algorithms, which are related to the individual classification criteria, CYCL and α-TUB were consistently ranked in the bottom five positions.
The consensus ranking obtained with RefFinder (Table 5), classified MDH as the most stable gene for all comparison groups, therefore constituting a good choice for the studies of Coffea spp. under future estimated conditions of increased [CO 2 ] and temperature. Furthermore, ACT and S15 were always ranked in the top five positions, while α-TUB and CYCL were amongst the least stable genes, as previously found in previous tests (Tables 3, 4). MDH and ACT have been previously described as being amongst the reference genes with higher transcript stability for RT-qPCR studies in Coffea spp. subjected to N starvation, salt and heat stresses (Carvalho et al., 2013) and to cold and drought stresses (Goulao et al., 2012), respectively. However, there was no consensus for the majority of candidate reference genes, i.e., the expression stability strongly varied according to  Table 1). Boxes indicate the 25th and 75th percentiles; lines across the box represent the median Ct-values; whisker caps represent the minimum and maximum values; and spots represent the outliers. the environmental stressful condition, again highlighting the importance of a proper and accurate selection of reference genes for each single experimental condition. For all comparison groups at V2/V3 the Pairwise variation (V) value was already below 0.15 (Figure 2), indicating that the use of two genes guarantee a strong accuracy for the normalization of RT-qPCR data. It is noteworthy to mention that for all treatments the pairwise variation values were <0.10 (except for the use of just one reference gene), indicating low expression instability of the selected genes tested for these environmental conditions.
In order to validate the consensus ranking, we have examined the relative expression of RLS, encoding the large-subunit of RuBisCO. The choice of this enzyme was related to its central role in the C-assimilation pathway, considering that improved crop photosynthetic efficiency under high temperature might be achieved through the replacement of a number of candidate amino acid substitutions to improve RuBisCO performance (Orr et al., 2016). Moreover, the RLS gene was selected based on the fact that, according to structural-functional analyses, most of the important amino acid residues necessary for catalysis are present in the large subunits (Morita et al., 2014). The consensus ranking established by RefFinder identified MDH, UBQ2, and DNAJ as the best stable genes to temperature, whereas a-TUB, CYCL and EF-1A where the least stable ones. Given that the increase in atmospheric temperature has proven to be the most disturbing variable to metabolism performance, mineral nutrition and acclimation limitation under predicted future environmental conditions (Martins et al., 2014Rodrigues et al., 2016), with particular strong negative implications to RuBisCO catalytic ability , we decided to use the most/least stable genes to temperature for RLS validation.
The expression profiles were normalized with the two and three most and least stable genes (Tables 6, 7), being observed that the expression patterns of RLS varied with the comparison For sake of simplicity, only significant differences were indicated (***P < 0.001; **P < 0.01; *P < 0.05) related to the BestKeeper test.
Frontiers in Plant Science | www.frontiersin.org FIGURE 2 | Pairwise variation analysis to determine the optimal number of reference genes for RT-qPCR data normalization.
group. Noticeably, due to the high SD of Ct values in many cases variations in RLS expression along the stress imposition were not reported as statistically significant. The patterns of RLS expression were nearly the same in RT-qPCR reactions calibrated with two or three genes. Independently of the stability of RLS mRNA, the increase in temperature generally resulted in upregulation of RLS expression at standard [CO 2 ] of 380 µL CO 2 L −1 in all genotypes, except in IPR108 at 31/25 • and 37/30 • C where RLS expression was nearly the same as the control. For the 700 µL CO 2 L −1 comparison group, the use of the most stable reference genes produced different results from those obtained with the least stable genes at 25/20 • and 31/25 • C. Thus, when reactions were normalized with the most stable genes, a consistent decrease in the steady-state amounts of RLS mRNA (which is the net result of transcription and degradation) was observed at these temperatures in all genotypes, while the use of the least stable genes indicated a positive regulation of RLS expression. For 37/30 • C and 42/34 • C, normalization of RT-qPCR with the most or the least stable genes produced similar up-regulated expression patterns of RLS. Notably, such an upregulation occurred when RuBisCO activity begin to decrease (37/30 • C) or was severely depressed (42/34 • C) . Such apparent discrepancy was also noted in chloroplastic ascorbate peroxidase (APX). APX activity was the most reduced (to about 10% of the control) at 42/34 • C among the studied enzymes, contrasting to the largest increase in APX transcripts amongst all studied genes . Such large increases in gene transcripts may be related to an attempt of the plants to compensate the severe decline of enzyme activity, and further highlights the need of integrated experiments that include gene expression and enzyme activity studies. Given the facts described above, the reference genes used in this study were quite stable under the applied experimental conditions, as pairwise variation values were <0.10. The results obtained in the groups 700-25/20 • C and 700-31/25 • C call up for the importance to select the best reference genes. Altogether, despite the low instability of a-TUB, CYCL, and EF-1A, the use of reference genes showing a high stability of their mRNA amount is the best option for the correct interpretation of data from transcriptional studies (Roy et al., 2011;Kozera and Rapacz, 2013). Also importantly, this corroborates with previous studies in Coffea sp., which indicate that the selection of reference genes For sake of simplicity, only significant differences were indicated (***P < 0.001; **P < 0.01; *P < 0.05).

CONCLUSIONS
In summary, our study reinforces the importance of the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines (Bustin et al., 2009) to normalize RT-qPCR data. The consensus ranking obtained with RefFinder showed that the transcript stability of the evaluated reference candidate genes changed according to the conditions/interactions, although MDH was the best one in all treatments, therefore, constituting a good common basis for RT-qPCR data normalization. However, the Pairwise variation analysis showed the need to use two reference genes, with the second most stable gene changing according to the genotypes and imposed environmental conditions: S15 for genotype and [CO 2 ], UBQ2 for temperature, EJF-4A for multiple stress, and ACT for multiple stress. Using the expression profiles of coffee RLS, results from the in silico aggregation and experimental validation For sake of simplicity, only significant differences were indicated (***P < 0.001; **P < 0.01; *P < 0.05).
confirmed that the use of two reference genes is an adequate procedure to normalize RT-qPCR data.

AUTHOR CONTRIBUTIONS
According to their competences, all authors contributed transversally to the several stages of the work, including its design, data acquisition, analysis and interpretation, critically review of the manuscript, and approval of the submitted version. Furthermore, they agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work were appropriately investigated and resolved.