Impact Factor 4.298

The 1st most cited journal in Plant Sciences

Original Research ARTICLE

Front. Plant Sci., 07 March 2017 | https://doi.org/10.3389/fpls.2017.00307

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

  • 1Plant-Environment Interactions and Biodiversity Lab (PlantStress&Biodiversity), Linking Landscape, Environment, Agriculture and Food, Departmento de Recursos Naturais, Ambiente e Território, Instituto Superior de Agronomia, Universidade de Lisboa (ULisboa), Oeiras, Portugal
  • 2Programa de Pós-Graduação em Genética e Melhoramento, Centro de Ciências Agrárias e Engenharias, Universidade Federal do Espírito Santo, Alegre, Brazil
  • 3Setor Fisiologia Vegetal, Centro de Ciências e Tecnologias Agropecuárias, Universidade Estadual Norte Fluminense-Darcy Ribeiro, Rio de Janeiro, Brazil
  • 4Departmento de Ciências Agrárias e Biológicas, Centro Universitário Norte do Espírito Santo, Universidade Federal Espírito Santo, São Mateus, Brazil
  • 5GeoBioTec, Departmento de Ciências da Terra, Faculdade de Ciências e Tecnologia, Universidade NOVA de Lisboa, Monte da Caparica, Portugal
  • 6Departmento de Biologia Vegetal, Universidade Federal Viçosa, Viçosa, Brazil

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.

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 [CO2] and enhanced temperatures predicted along this century. In fact, depending on anthropogenic greenhouse gas emission scenarios, [CO2] might reach between 421 and 936 μL CO2 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, 2014). However, it was recently shown that increased [CO2] in the atmosphere strengthens the coffee plant by reinforcing its photosynthetic performance (Ramalho et al., 2013; Ghini et al., 2015; DaMatta et al., 2016). Furthermore, enhanced [CO2] 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 (Rodrigues et al., 2016), should be reconsidered. In addition, results from Rodrigues et al. (2016) strongly suggest the need to study the single and superimposed effects of elevated [CO2] 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://coffee-genome.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 (Wang et al., 2016). In addition, the species-dependent 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 [CO2] and heat were considered.

Materials and Methods

Plant Material and Growth Conditions

Plant materials and experimental design was previously described in detail (Rodrigues et al., 2016). 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 [CO2]) into walk-in growth chambers (EHHF 10,000, ARALAB, Portugal) differing in air [CO2] supply: 380 μL CO2 L−1 (380-plants) or 700 μL CO2 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., 2014, 2016; Rodrigues 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 [CO2] (380, 700 μL CO2 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® DNA-free™ 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.

TABLE 1
www.frontiersin.org

Table 1. Primer sequences and amplicon characteristics for each of the 10 candidate reference genes under evaluation.

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).

Selection of Reference Gene Sequences and Primer Design

A set of 10 candidate genes were selected, comprising several frequently used reference genes, based on previous reports in (a) coffee plants (Cruz et al., 2009; Goulao et al., 2012; Carvalho et al., 2013): actin (ACT), glyceraldehyde 3-phosphate dehydrogenase (GAPDH), eukaryotic initiation factor 4α (Elf4A), elongation factor 1α (EF1A), cyclophilin (CYCL), and malate dehydrogenase (MDH); and (b) other plants (Cseke et al., 2009; El Kelish et al., 2014; Watanabe et al., 2014; Yan et al., 2014): 40S ribosomal protein subunit 15A (S15), DNAJ-like (DNAJ), ubiquitin-conjugating enzyme E2(UBQ2), and a alpha tubulin (α-TUB). The corresponding gene sequences were retrieved from published literature and the NCBI databases (Cruz et al., 2009; Vidal et al., 2010; Mondego et al., 2011; Goulao et al., 2012; Carvalho et al., 2013). Primers were designed using Primer3 software (Rozen and Skaletsky, 2000). The length of the primers was set to be between 19 and 26 bp, with a GC content ranging from 45 to 60% and a melting temperature (Tm-value) between 62° and 65°C. Amplicon length ranges were set to be between 100 and 150 bp. The probability of formation of hairpin structures and primer dimerization was subsequently checked using the Oligo Calculator (version. 3.26) algorithm (Kibbe, 2007). Primer sequences are given in Table 1.

Quantitative Real-Time Polymerase Chain Reaction Conditions

RT-qPCR reactions were performed in 96-well plates using iQ™ SYBR® 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™ SYBR® Green Supermix. Reactions were run on the iQ™ 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: E (%) = (10−(1/slope)−1) × 100 (da 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 web-based 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 (Vn/Vn+1) of two consecutively serial log2-transformed normalization factor (NF) ratios, based on stepwise addition of the subsequent more stable reference gene (NFn relative to NFn+1 reference genes; Vandesompele et al., 2002). The SD of the array generated by the log2-transformed NF ratios was calculated for each NF combination (Vn/Vn+1). Vn/Vn+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 [CO2] (380, 700 μL CO2 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, [CO2] and temperature, multiple stress interaction ([CO2] and temperature) and total stress interaction (genotype, [CO2] 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 (R2) 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).

TABLE 2
www.frontiersin.org

Table 2. Efficiency values and correlation coefficient for each candidate reference genes under evaluation.

FIGURE 1
www.frontiersin.org

Figure 1. Box plot representing the expression profiling variations of each candidate reference gene. The primer pairs of each gene were used to examine all samples (n = 24; 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.

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 [CO2], 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 [CO2] and multiple stress. DNAJ and TUB were the less stable genes for [CO2] (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, [CO2] 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 [CO2] 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).

TABLE 3
www.frontiersin.org

Table 3. Ranking of candidate reference genes according to its coefficient of variation (CV%), considering the variables genotype, temperature, [CO2] and their interaction.

TABLE 4
www.frontiersin.org

Table 4. Stability of candidate reference genes according to GeNorm, NormFinder and BestKeeper, considering the variables genotype, temperature, [CO2], and their interaction.

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 [CO2] 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 [CO2] 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 the environmental stressful condition, again highlighting the importance of a proper and accurate selection of reference genes for each single experimental condition.

TABLE 5
www.frontiersin.org

Table 5. Overall ranking of the most stable genes within each treatment made through RefFinder program, considering the variables genotype, temperature, [CO2] and their interaction.

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.

FIGURE 2
www.frontiersin.org

Figure 2. Pairwise variation analysis to determine the optimal number of reference genes for RT-qPCR data normalization.

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., 2014, 2016; Rodrigues et al., 2016), with particular strong negative implications to RuBisCO catalytic ability (Rodrigues et al., 2016), 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 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 up-regulation of RLS expression at standard [CO2] of 380 μL CO2 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 CO2 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 up-regulation occurred when RuBisCO activity begin to decrease (37/30°C) or was severely depressed (42/34°C) (Rodrigues et al., 2016). 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 (Martins et al., 2016). 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.

TABLE 6
www.frontiersin.org

Table 6. Analysis of RLS expression by RT-qPCR, calibrated with the two most (MDH and UBQ2) and the two least (a-TUB and CYCL) genes for the temperature group/treatment.

TABLE 7
www.frontiersin.org

Table 7. Analysis of RLS expression by RT-qPCR, calibrated with the three most (MDH, UBQ2 and DNAJ) and the three least (a-TUB, CYCL and EF-1A) genes for the temperature group/treatment.

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 an accurate RT-qPCR normalization should be performed accordingly to each for specific condition (Cruz et al., 2009; Goulao et al., 2012; Carvalho et al., 2013; Yuyama et al., 2016).

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 [CO2], 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 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.

Funding

This work was supported by by national funds from Fundação para a Ciência e a Tecnologia through the projects PTDC/AGR-PRO/3386/2012, the research units UID/AGR/04129/2013 (LEAF) and UID/GEO/04035/2013 (GeoBioTec), as well through the grant SFRH/BPD/47563/2008 (AF) co-financed through the POPH program subsidized by the European Social Fund. Brazilian funding from CAPES (grants PDSE: 000427/2014-04, WR; 0343/2014-05, MM), CNPq and Fapemig (fellowships to FD, FP, and EC) are also greatly acknowledged.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank DELTA Cafés (Portugal) for technical assistance, and Profs. T. Sera (IAPAR) for supplying IPR108 material. Thanks are also due to M. Costa, and N. Duro, for technical help.

Abbreviations

RT-qPCR, quantitative real-time polymerase chain reaction.

References

Andersen, C. L., Jensen, J. L., and Orntoft, T. F. (2004). Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64, 5245–5250. doi: 10.1158/0008-5472.CAN-04-0496

PubMed Abstract | CrossRef Full Text | Google Scholar

Artico, S., Nardeli, S. M., Brilhante, O., Grossi-de-Sa, M. F., and Alves-Ferreira, M. (2010). Identification and evaluation of new reference genes in Gossypium hirsutum for accurate normalization of real time quantitative RT-PCR data. BMC Plant Biol. 10:49. doi: 10.1186/1471-2229-10-49

PubMed Abstract | CrossRef Full Text | Google Scholar

Bunn, C., Läderach, P., Ovalle-Rivera, O., and Kirschke, D. (2015). A bitter cup: climate change profile of global production of Arabica and Robusta coffee. Clim. Change 129, 89–101. doi: 10.1007/s10584-014-1306-x

CrossRef Full Text | Google Scholar

Bustin, S. A., Benes, V., Garson, J. A., Hellemans, J., Huggett, J., Kubista, M., et al. (2009). The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55, 611–622. doi: 10.1373/clinchem.2008.112797

PubMed Abstract | CrossRef Full Text | Google Scholar

Carvalho, K., Bespalhok Filho, J. C., dos Santos, T. B., de Souza, S. G., Vieira, L. G., Pereira, L. F., et al. (2013). Nitrogen starvation, salt and heat stress in coffee (Coffea arabica L.): identification and validation of new genes for qPCR normalization. Mol. Biotechnol. 53, 315–325. doi: 10.1007/s12033-012-9529-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cruz, F., Kalaoun, S., Nobile, P., Colombo, C., Almeida, J., Barros, L., et al. (2009). Evaluation of coffee reference genes for relative expression studies by quantitative real-time RT–PCR. Mol. Breed. 23, 607–616. doi: 10.1007/s11032-009-9259-x

CrossRef Full Text | Google Scholar

Cseke, L. J., Tsai, C.-J., Rogers, A., Nelsen, M. P., White, H. L., Karnosky, D. F., et al. (2009). Transcriptomic comparison in the leaves of two aspen genotypes having similar carbon assimilation rates but different partitioning patterns under elevated [CO2]. New Phytol. 182, 891–911. doi: 10.1111/j.1469-8137.2009.02812.x

CrossRef Full Text | Google Scholar

da Costa, M., Duro, N., Batista-Santos, P., Ramalho, J. C., and Ribeiro-Barros, A. I. (2015). Validation of candidate reference genes for qRT-PCR studies in symbiotic and non-symbiotic Casuarina glauca Sieb. ex Spreng. under salinity conditions. Symbiosis 66, 21–35. doi: 10.1007/s13199-015-0330-6

CrossRef Full Text | Google Scholar

DaMatta, F. M., Godoy, A. G., Menezes-Silva, P. E., Martins, S. C. V., Sanglard, L. M. V. P., Morais, L. E., et al. (2016). Sustained enhancement of photosynthesis in coffee trees grown under free-air CO2 enrichment conditions: disentangling the contributions of stomatal, mesophyll, and biochemical limitations. J. Exp. Bot. 67, 341–352. doi: 10.1093/jxb/erv463

PubMed Abstract | CrossRef Full Text | Google Scholar

DaMatta, F. M., and Ramalho, J. C. (2006). Impacts of drought and temperature stress on coffee physiology and production: a review. Theor. Exp. Plant Physiol. 18, 55–81. doi: 10.1590/S1677-04202006000100006

CrossRef Full Text | Google Scholar

Die, J. V., Román, B., Nadal, S., and González-Verdejo, C. I. (2010). Evaluation of candidate reference genes for expression studies in Pisum sativum under different experimental conditions. Planta 232, 145–153. doi: 10.1007/s00425-010-1158-1

PubMed Abstract | CrossRef Full Text | Google Scholar

El Kelish, A., Zhao, F., Heller, W., Durner, J., Winkler, J. B., Behrend, H., et al. (2014). Ragweed (Ambrosia artemisiifolia) pollen allergenicity: SuperSAGE transcriptomic analysis upon elevated CO2 and drought stress. BMC Plant Biol. 14:176. doi: 10.1186/1471-2229-14-176

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghini, R., Torre-Neto, A., Dentzien, A. F. M., Guerreiro-Filho, O., Lost, R., Patrício, F. R. A., et al. (2015). Coffee growth, pest and yield responses to free-air CO2 enrichment. Clim. Chang. 132, 307–320. doi: 10.1007/s10584-015-1422-2

CrossRef Full Text | Google Scholar

Goulao, L. F., Fortunato, A. S., and Ramalho, J. C. (2012). Selection of reference genes for normalizing quantitative real-time PCR gene expression data with multiple variables in Coffea spp. Plant Mol. Biol. Rep. 30, 741–759. doi: 10.1007/s11105-011-0382-6

CrossRef Full Text | Google Scholar

Gutierrez, L., Mauriat, M., Guenin, S., Pelloux, J., Lefebvre, J.-F., Louvet, R., et al. (2008b). The lack of a systematic validation of reference genes: a serious pitfall undervalued in reverse transcription-polymerase chain reaction (RT-PCR) analysis in plants. Plant Biotechnol. J. 6, 609–618. doi: 10.1111/j.1467-7652.2008.00346.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gutierrez, L., Mauriat, M., Pelloux, J., Bellini, C., and van Wuytswinkel, O. (2008a). Towards a systematic validation of references in real-time RT-PCR. Plant Cell 20, 1734–1735. doi: 10.1105/tpc.108.059774

PubMed Abstract | CrossRef Full Text | Google Scholar

ICO (2014a). Trade Statistics. Available online at: http://www.ico.org/trade_statistics.asp (Accessed July 2016)

ICO (2014b). World coffee trade (1963–2013): A Review of the Markets, Challenges and Opportunities Facing the Sector. ICC (International Coffee Council), 111–115 Rev. 1, 29p. Available online at: http://www.ico.org/show_news.asp?id=361 (Accessed July 2016)

Imai, T., Ubi, B. E., Saito, T., and Moriguchi, T. (2014). Evaluation of reference genes for accurate normalization of gene expression for real time-quantitative PCR in Pyrus pyrifolia using different tissue samples and seasonal conditions. PLoS ONE 9:e86492. doi: 10.1371/journal.pone.0086492

PubMed Abstract | CrossRef Full Text | Google Scholar

IPCC (2013). “Climate change 2013,” in The Physical Science Basis. Summary for Policymakers, Technical Summary and Frequent Asked Questions. Part of the Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, eds T.F. Stocker, D. Qin, G. K. Plattner, M. M. B. Tignor, S. K. Allen, and J. Boschung, et al. (Intergovernmental Panel on Climate Change), 203.

IPCC (2014). “Climate change 2014,” in Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, eds O. Edenhofer, R. Pichs-Madruga, Y. Sokona, E. Farahani, S. Kadner, and K. Seyboth, et al. (Cambridge, New York, NY: Cambridge University Press), 1435.

Kibbe, W. A. (2007). OligoCalc: an online oligonucleotide properties calculator. Nucl. Acids Res. 35, W43–W46. doi: 10.1093/nar/gkm234

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozera, B., and Rapacz, M. (2013). Reference genes in real-time PCR. J. Appl. Gen. 54, 391–406. doi: 10.1007/s13353-013-0173-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M. Y., Wang, F., Jiang, Q., Wang, G. L., Tian, C., and Xiong, A. S. (2016). Validation and comparison of reference genes for qPCR normalization of celery (Apium graveolens) at different development stages. Front. Plant Sci. 7:313. doi: 10.3389/fpls.2016.00313

PubMed Abstract | CrossRef Full Text | Google Scholar

Llanos, A., François, J. M., and Parrou, J. L. (2015). Tracking the best reference genes for RT-qPCR data normalization in filamentous fungi. BMC Genomics. 16:71. doi: 10.1186/s12864-015-1224-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Magrach, A., and Ghazoul, J. (2015). Climate and pest-driven geographic shifts in global coffee production: implications for forest cover, biodiversity and carbon storage. PLoS ONE 10:e0133071. doi: 10.1371/journal.pone.0133071

PubMed Abstract | CrossRef Full Text | Google Scholar

Mallona, I., Lischewski, S., Weiss, J., Hause, B., and Egea-Cortines, M. (2010). Validation of reference genes for quantitative realtime PCR during leaf and flower development in Petunia hybrida. BMC Plant Biol. 10:4. doi: 10.1186/1471-2229-10-4

CrossRef Full Text | Google Scholar

Martins, L. D., Tomaz, M. A., Lidon, F. C., DaMatta, F. M., and Ramalho, J. C. (2014). Combined effects of elevated [CO2] and high temperature on leaf mineral balance in Coffea spp. plants. Clim. Chang. 126, 365–379. doi: 10.1007/s10584-014-1236-7

CrossRef Full Text | Google Scholar

Martins, M. Q., Rodrigues, W. P., Fortunato, A. S., Leitão, A. E., Rodrigues, A. P., Pais, I. P., et al. (2016). Protective response mechanisms to heat stress in interaction with high [CO2] conditions in Coffea spp. Front. Plant Sci. 7:947. doi: 10.3389/fpls.2016.00947

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondego, J. M., Vidal, R. O., Carazzolle, M. F., Tokuda, E. K., Parizzi, L. P., Costa, G. G., et al. (2011). An EST-based analysis identifies new genes and reveals distinctive gene expression features of Coffea arabica and Coffea canephora. BMC Plant Biol. 11:30. doi: 10.1186/1471-2229-11-30

PubMed Abstract | CrossRef Full Text | Google Scholar

Morita, K., Hatanaka, T., Misoo, S., and Fukayama, H. (2014). Unusual small subunit that is not expressed in photosynthetic cells alters the catalytic properties of rubisco in rice. Plant Physiol. 164, 69–79. doi: 10.1104/pp.113.228015

CrossRef Full Text | Google Scholar

Niu, X., Qi, J., Zhang, G., Xu, J., Tao, A., Fang, P., et al. (2015). Selection of reliable reference genes for quantitative real-time PCR gene expression analysis in Jute (Corchorus capsularis) under stress treatments. Front. Plant Sci. 6:848. doi: 10.3389/fpls.2015.00848

PubMed Abstract | CrossRef Full Text | Google Scholar

Orr, D. J., Alcβntara, A., Kapralov, M. V., Andralojc, P. J., Carmo-Silva, E., and Parry, M. A. J. (2016). Surveying rubisco diversity and temperature response to improve crop photosynthetic efficiency. Plant Physiol. 172, 707–717. doi: 10.1104/pp.16.00750

PubMed Abstract | CrossRef Full Text | Google Scholar

Petriccione, M., Mastrobuoni, F., Zampella, L., and Scortichini, M. (2015). Reference gene selection for normalization of RT-qPCR gene expression data from Actinidia deliciosa leaves infected with Pseudomonas syringae pv. actinidiae. Sci. Rep. 5:16961. doi: 10.1038/srep16961

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfaffl, M. W., Tichopad, A., Prgomet, C., and Neuvians, T. P. (2004). Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestkeeper - excel-based tool using pair-wise correlations. Biotechnol. Lett. 26, 509–515. doi: 10.1023/B:BILE.0000019559.84305.47

PubMed Abstract | CrossRef Full Text | Google Scholar

Poncet, V., Rondeau, M., Tranchant, C., Cayrel, A., Hamon, S., De Kochko, A., et al. (2006). SSR mining in coffee tree EST databases: potential use of EST-SSRs as markers for the Coffea genus. Mol. Genet. Genomics 276, 436–449. doi: 10.1007/s00438-006-0153-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramalho, J. C., Rodrigues, A. P., Semedo, J. N., Pais, I. P., Martins, L. D., Simões-Costa, M. C., et al. (2013). Sustained photosynthetic performance of Coffea spp. under long-term enhanced [CO2]. PLoS ONE 8:e82712. doi: 10.1371/journal.pone.0082712

CrossRef Full Text | Google Scholar

R Core Team (2016). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. Available online at: https://www.R-project.org (Accessed November 2016).

Rodrigues, W. P., Martins, M. Q., Fortunato, A. S., Rodrigues, A. P., Semedo, J. N., Simões-Costa, M. C., et al. (2016). Long-term elevated air [CO2] strengthens photosynthetic functioning and mitigates the impact of supra-optimal temperatures in tropical Coffea arabica and C. canephora species. Glob. Change Biol. 22, 415–431. doi: 10.1111/gcb.13088

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, S. J., Tucker, E. J., and Tester, M. (2011). Genetic analysis of abiotic stress tolerance in crops. Curr. Opin. Plant Biol. 14, 232–239. doi: 10.1016/j.pbi.2011.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Rozen, S., and Skaletsky, H. (2000). “Bioinformatics methods and protocols: methods,” in Molecular Biology, Chapter 20, Primer3 on the WWW for general users and for biologist programmers, eds S. Krawetz and S. Misener (Totowa, NJ: Humana), 365–386.

Tian, C., Jiang, Q., Wang, F., Wang, G. L., Xu, Z. S., and Xiong, A. S. (2015). Selection of suitable reference genes for qPCR normalization under abiotic stresses and hormone stimuli in carrot leaves. PLoS ONE 10:e0117569. doi: 10.1371/journal.pone.0117569

PubMed Abstract | CrossRef Full Text | Google Scholar

Vandesompele, J., De Preter, K., Pattyn, F., Poppe, B., Van Roy, N., De Paepe, A., et al. (2002). Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3, research0034.1–research0034.11. doi: 10.1186/gb-2002-3-7-research0034

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidal, R. O., Mondego, J. M. C., Pot, D., Ambrósio, A. B., Andrade, A. C., Pereira, L. F. P., et al. (2010). A high-throughput data mining of single nucleotide polymorphisms in coffea species expressed sequence tags suggests differential homeologous gene expression in the allotetraploid Coffea arabica. Plant Physiol. 154, 1053–1066. doi: 10.1104/pp.110.162438

PubMed Abstract | CrossRef Full Text | Google Scholar

Vieira, L. G. E., Andrade, A. C., Colombo, C. A., Moraes, A. A. H., Metha, A., Oliveira, A. C., et al. (2006). Brazilian coffee genome project: an EST-based genomic resource. Theor. Exp. Plant Physiol. 18, 95–108. doi: 10.1590/S1677-04202006000100008

CrossRef Full Text | Google Scholar

Waller, J. M., Bigger, M., and Hillocks, R. J. (2007). “Coffee pests, diseases and their management,” in World Coffee Production, Chapter 2. eds J. M. Waller, M. Bigger, and R.J. Hillocks (Egham; Surrey, BC: CAB International,) 17–40.

Google Scholar

Wang, G. L., Tian, C., Jiang, Q., Xu, Z. S., Wang, F., and Xiong, A. S. (2016). Comparison of nine reference genes for real-time quantitative PCR in roots and leaves during five developmental stages in carrot (Daucus carota L.). J. Hortic. Sci. Biotechnol. 91, 264–270. doi: 10.1080/14620316.2016.1148372

CrossRef Full Text | Google Scholar

Watanabe, C. K., Sato, S., Yanagisawa, S., Uesono, Y., Terashima, I., and Noguchi, K. (2014). Effects of elevated CO2 on levels of primary metabolites and transcripts of genes encoding respiratory enzymes and their diurnal patterns in Arabidopsis thaliana: possible relationships with respiratory rates. Plant Cell Physiol. 55, 341–357. doi: 10.1093/pcp/pct185

PubMed Abstract | CrossRef Full Text | Google Scholar

Willems, E., Leyns, L., and Vandesompele, J. (2008). Standardization of real-time PCR gene expression data from independent biological replicates. Anal. Biochem. 379, 127–129. doi: 10.1016/j.ab.2008.04.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, F., Xiao, P., Chen, D., Xu, L., and Zhang, B. (2012). miRDeepFinder: a miRNA analysis tool for deep sequencing of plant small RNAs. Plant Mol. Biol. 80, 75–84. doi: 10.1007/s11103-012-9885-2

CrossRef Full Text | Google Scholar

Yan, X., Dong, X., Zhang, W., Yin, H., Xiao, H., Chen, P., et al. (2014). Reference gene selection for quantitative Real-Time PCR normalization in Reaumuria soongorica. PLoS ONE 9:e104124. doi: 10.1371/journal.pone.0104124

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Pan, H., Noland, J. E., Zhang, D., Zhang, Z., Liu, Y., et al. (2015). Selection of reference genes for RT-qPCR analysis in a predatory biological control agent, Coleomegilla maculata (Coleoptera: Coccinellidae). Sci. Rep. 5:18201. doi: 10.1038/srep18201

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, X., Zhang, F., Tao, Y., Song, S., and Fang, J. (2015). Reference gene selection for quantitative real-time PCR normalization in different cherry genotypes, developmental stages and organs. Sci. Hortic. 181, 182–188. doi: 10.1016/j.scienta.2014.10.027

CrossRef Full Text | Google Scholar

Yuyama, P. M., Reis Júnior, O., Ivamoto, S. T., Domingues, D. S., Carazzolle, M. F., Pereira, G. A., et al. (2016). Transcriptome analysis in Coffea eugenioides, an Arabica coffee ancestor, reveals differentially expressed genes in leaves and fruits. Mol. Genet. Genomics 291, 323–336. doi: 10.1007/s00438-015-1111-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: climate changes, coffee, increased air [CO2], global warming, heat stress, normalization of transcriptomic studies, quantitative real-time PCR, reference genes

Citation: Martins MQ, Fortunato AS, Rodrigues WP, Partelli FL, Campostrini E, Lidon FC, DaMatta FM, Ramalho JC and Ribeiro-Barros AI (2017) 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. Front. Plant Sci. 8:307. doi: 10.3389/fpls.2017.00307

Received: 22 December 2016; Accepted: 20 February 2017;
Published: 07 March 2017.

Edited by:

Bernd Mueller-Roeber, University of Potsdam, Germany

Reviewed by:

Andrew Wood, Southern Illinois University Carbondale, USA
Dirk K. Hincha, Max Planck Institute of Molecular Plant Physiology (MPG), Germany

Copyright © 2017 Martins, Fortunato, Rodrigues, Partelli, Campostrini, Lidon, DaMatta, Ramalho and Ribeiro-Barros. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: José C. Ramalho, cochichor@mail.telepac.pt; cochichor@isa.ulisboa.pt

These authors have contributed equally to the work.