The Greater Phenotypic Homeostasis of the Allopolyploid Coffea arabica Improved the Transcriptional Homeostasis Over that of Both Diploid Parents

Abstract

Polyploidy impacts the diversity of plant species, giving rise to novel phenotypes and leading to ecological diversification. In order to observe adaptive and evolutionary capacities of polyploids, we compared the growth, primary metabolism and transcriptomic expression level in the leaves of the newly formed allotetraploid Coffea arabica species compared with its two diploid parental species ( Coffea eugenioides and Coffea canephora ), exposed to four thermal regimes (TRs; 18–14, 23–19, 28–24 and 33–29°C). The growth rate of the allopolyploid C. arabica was similar to that of C. canephora under the hottest TR and that of C. eugenioides under the coldest TR. For metabolite contents measured at the hottest TR, the allopolyploid showed similar behavior to C. canephora , the parent which tolerates higher growth temperatures in the natural environment. However, at the coldest TR, the allopolyploid displayed higher sucrose, raffinose and ABA contents than those of its two parents and similar linolenic acid leaf composition and Chl content to those of C. eugenioides . At the gene expression level, few differences between the allopolyploid and its parents were observed for studied genes linked to photosynthesis, respiration and the circadian clock, whereas genes linked to redox activity showed a greater capacity of the allopolyploid for homeostasis. Finally, we found that the overall transcriptional response to TRs of the allopolyploid was more homeostatic compared with its parents. This better transcriptional homeostasis of the allopolyploid C. arabica afforded a greater phenotypic homeostasis when faced with environments that are unsuited to the diploid parental species.

Introduction

Allopolyploidization has long been recognized as a major mechanism of adaptation and speciation in plants. Allopolyploid species, defined as interspecific hybrids in which chromosome sets are doubled, lead to permanent fixation of heterozygosity and hybrid vigor ( Chen 2010 ). Fixed heterozygosity is well documented as a source of genetic variation and phenotypic variability in polyploids ( Feldman and Levy 2005 ) and could contribute to a greater ecological tolerance of polyploids compared with their parents. Many phenotypic changes with a high adaptive value, concerning vigor, biomass and adaptation to biotic and abiotic stress, would enable the polyploids to occupy new ecological niches and seem to be conducive to their use in agriculture ( Osborn et al. 2003 , Chen 2007 ). A study by Ni et al. (2009) focusing on allopolyploids of Arabidopsis sp. revealed some better performances in allopolyploid individuals which could be assimilated to the heterosis seen in hybrids, but with the particularity of being fixed in the allopolyploidy context. In addition, allopolyploidy is often associated with enhanced stress tolerance. For example, the Spartina anglica allopolyploid is adapted to varied environments and also displayed better tolerance than its parental species to reducing conditions and to sulfite-rich sediments ( Ainouche et al. 2004 , Maricle et al. 2006 , Ainouche et al. 2009 ).

In order to observe adaptive and evolutionary capacities in allopolyploids, we used Coffea arabica , a tropical perennial species, as our study model. Coffea arabica is an allopolyploid species (2 n = 4 x = 44) resulting from the hybridization between Coffea eugenioides and Coffea canephora or ecotypes related to those diploid species ( Lashermes et al. 1999 ). The allopolyploid speciation of C. arabica took place in relatively recent times. The analysis of nucleotide substitutions dates the origin of the allotetraploid species C. arabica in a time interval spanning 10,000 or 150,000 years BP ( Cenci et al. 2012 ). In the wild state, the territory of C. canephora stretches from Central Africa to the Gulf of Guinea and Uganda. Compared with the vast territory occupied by C. canephora , the territories of C. arabica and C. eugenioides are more limited. The territory of wild C. eugenioides stretches from east to south-east Africa (Uganda, Tanzania, Kenya, Rwanda and Burundi). For its part, C. arabica comes from a limited region of south-western Ethiopia and the Boma Plateau of Sudan where no other Coffea species is represented. Thus, these three Coffea species also grow at different elevations and temperatures. The wild C. arabica species is found in high altitude tropical forests, between 1,200 and 1,950 m ( Davis et al. 2006 ), with a mean annual temperature optimum of 18–21°C ( Damatta and Ramalho 2006 ). The wild C. canephora species comes from low to medium altitude equatorial rainforests, between 250 and 1,500 m ( Davis et al. 2006 ), with a mean annual temperature optimum of 22–26°C ( Damatta and Ramalho 2006 ). Lastly, the wild C. eugenioides species is found in high altitude forests between 1,000 and 2,000 m, with a mean annual temperature optimum of 18–23°C ( Davis et al. 2006 ).

Temperature provides a good framework for testing homeostasis (vs. phenotypic plasticity), because many traits are highly sensitive to changes in thermal conditions ( Bita and Gerats 2013 ). In this study, we tried to establish performances of C. arabica compared with its parental diploid species under thermal regimes (TRs), representing a possible range of intertropical TRs. We compared the growth, primary metabolism and patterns of transcriptomic expression of these three species exposed to different TRs. Following the definition of Schlitting and Smith (2002) and Nicotra et al. (2010) , phenotypic plasticity is construed to encompass a diversity of phenomena spanning several hierarchical levels of organization and assumed to be mediated at the molecular level. Assuming that homeostasis (vs. phenotypic plasticity) is mediated at the molecular level ( Schlichting and Smith 2002 , Nicotra et al. 2010 ), we measured the level of gene expression in leaves taking an RNA sequencing (RNAseq) approach, and an index that globally measured the plasticity of the transcriptomes of the three species depending on the environment is calculated. Gene expression phenotypes can be also measured as a reaction norm in response to environmental variations ( Hodgins-Davis and Townsend 2009 ). The genes were classed according to their stability (i.e. plasticity) depending on the thermal environment and/or the species. Groups of genes corresponding to photosynthesis, coding for pigments and proteins involved in electron transport or genes coding for related aspects of carbon metabolism were particularly observed. Indeed, according to Ni et al. (2009) , an up-regulation of these genes should be anticipated in the allopolyploid. The response of genes involved in the reactive oxygen species (ROS) scavenging network may also be of interest especially when growth temperatures become stressful for plants (18–14 and 33–29°C). Lastly, we particularly studied some key genes of the circadian clock, since alteration of the circadian rhythm would appear to play a major role in growth vigor and increased biomass in Arabidopsis allotetraploids ( Ni et al. 2009 ) and because they regulate the abundance of proteins involved in many primary metabolic pathways ( Hotta et al . 2007 ).

The purpose of this study was 3-fold. On the one hand, we studied a relatively recent natural allopolyploid in evolutionary terms, but with a sufficient number of generations to be stable, well regulated at the meiosis level and subjected to natural selection, whereas most studies were conducted with neosynthetic allopolyploid models (often poorly regulated at the meiosis level or/and presenting a high level of hybrid vigor). On the other hand, we studied a perennial tropical species, whereas annual temperate species were most often used in studies on phenotypic plasticity depending on TRs. Lastly, we used four contrasting tropical TRs, whereas most studies were carried out in a single environment.

Results

The allopolyploid growth rate was similar to those of its parental diploid species depending on thermal regimes

We compared the slopes of the thermal reaction norms among the allopolyploid species C. arabica and its diploid parental species ( C. canephora and C. eugenioides ) ( Supplementary Table S1 ). A relative growth rate was calculated from the difference between the dry weights estimated before and after the period in phytotrons (at T0 and T11 weeks). The growth of the three species differed depending on the TRs. Under the coldest TR (18–14°C), C. canephora growth appeared to be affected but no mortality was observed ( Fig. 1 A). Conversely, C. eugenioides growth was greatly affected by the hottest TR (33–29°C), with a mortality amounting to 52%. We considered that this TR was highly stressful for both species. On the other hand, the allopolyploid growth seemed stable whatever the TR, and no mortality was observed ( Fig. 1 A). Its growth rate was similar to that of C. canephora under the hottest TR (33–29°C) and to that of C. eugenioides under the coldest TR (18–14°C).

Fig. 1

Influence of temperature on growth and coffee leaf composition of three Coffea species. (A) Growth rate, (B) Chl, (C) total sugars, (D) sucrose, (E) stacchyose and (F) raffinose. Open circles, filled squares and open triangles are for Coffea canephora , Coffea arabica and Coffea eugenioides , respectively . For growth increase, the data are the means ± SEM ( n = 15). For Chl and sugar contents, values and error bars are means ± SEM of three experimental replicates. Different letters indicate significant differences (a,b and c; P < 0.05) between species for each thermal regime. ns, non-significant differences.

A dramatic drop of Coffea canephora Chl content at the coldest thermal regime (18–14°C)

Whatever the species considered, the Chl content declined under the coldest TR (18–14°C; Fig. 1 B). A dramatic drop of Chl content was observed for C. canephora under the coldest TR and some photobleaching symptoms were found on fully expanded leaves. The allopolyploid Chl content was similar to that of C. eugenioides under the two coldest TRs (18–14 and 23–19°C) and to that of C. canephora under the two hottest TRs (28–24 and 33–29°C).

Sucrose and raffinose family oligosaccharide contents: original performance of the Coffea arabica allopolyploid

Among the total sugar content of leaves ( Fig. 1 C) sampled at Zeitgeber time = 6 h (ZT6), 60–81% is identified as sucrose ( Fig. 1 D), 11–35% as glucose and fructose (data not shown), 5–7% as raffinose and stachyose ( Fig. 1 E, F) and <5% as other sugars such as maltose, inositol and sorbose (data not shown).

The allopolyploid sucrose content was stable for the three coldest TRs (18–14, 23–19 and 28–24°C; Fig. 1 D). That value fell by more than half under the hottest TR (33–29°C). Under the 23–19°C TR, there did not appear to be any differences between the three species. Under the 18–14 and 28–24°C TRs, the allopolyploid species displayed a higher sucrose content than those of its two parental species, which would seem to reflect more efficient photosynthesis or a signaling role. For the hottest TR (33–29°C), no difference was found between C. arabica and C. canephora. The raffinose family oligosaccharides (RFOs), represented in this study by raffinose and stachyose ( Fig. 1 E, F), are synthesized from sucrose. For stachyose content, there did not appear to be any significant differences between both parental diploid species, and content was relatively stable between the TRs ( Fig. 1 E). For the three coldest TRs (18–14, 23–19 and 28–24°C), the allopolyploid stachyose content was stable and significantly higher than those of both parental diploid species. At the hottest TR (33–29°C), the allopolyploid stachyose content was lower than for the other TRs and similar to the C. canephora stachyose content. The C. canephora raffinose ( Fig. 1 F) content was stable for the three hottest TRs (23–19, 28–24 and 33–29°C) and increased under the coldest TR (18–14°C). The tendency seemed to be reversed in C. eugenioides , since an increase in contents was found when the temperatures rose. The allopolyploid raffinose content was stable and higher than those of the two parental diploid species for the three coldest TRs (18–14, 23–19 and 28–24°C). Under the hottest TR (33–29°C), the C. arabica raffinose content dropped substantially and was similar to that of C. canephora .

The leaf linolenic acid content of the allopolyploid was similar to those of its parental diploid species depending on thermal regimes

In the three Coffea species studied, palmitic acid (C16:0) and linolenic acid (C18:3) are the two most abundant fatty acids. Variations in the relative abundance of these fatty acids were found for all species and TRs ( Fig. 2 A, B). Depending on the TRs, percentages ranging from 28% to 50% were found for C16: 0 and from 11% to 39% for C18:3. For all the TRs, the leaf C16:0 content was higher in C. eugenioides (42.5–50.2%) than in C. canephora (>30%), whereas the C18:3 content was greater in C. canephora (27–39.1%) than in C. eugenioides (11.4–29.8%). Under the coldest TR (18–14°C), the C18:3 percentage increased for C. canephora and C. eugenioides . For the allopolyploid C. arabica , there was a high variability in the C16:0 and C18:3 contents depending on TRs. Under the two hottest TRs (28–24 and 33–29°C), it was found that the C. arabica C16:0 and C18:3 percentages were similar to those of C. canephora . For the two coldest TRs (18–14 and 23–19°C), they were intermediate between the two parental diploid species for C16:0 and similar to those of C. eugenioides for C18:3. The other five fatty acids (palmitoleic acid, stearic acid, oleic acid, linoleic acid and arachidic acid) accounted for 25–39% of the total (data not shown). For these, no significant variations were found, either between species or between TRs.

Fig. 2

Influence of temperature on coffee leaf composition of three Coffea species. (A) Palmitic acid (C16:0), (B) linolenic acid (C18:3), (C) (D) ABA–glucose ester conjugate (ABAGE), (E) dihydrophaseic acid (DPA) and (F) ABA. Open circles, filled squares and open triangles are for Coffea canephora , Coffea arabica and Coffea eugenioides , respectively . Values and error bars are means ± SE of three experimental replicates. Different letters indicate significant differences (a,b and c; P < 0.05) between species for each thermal regime. ns, non-significant differences. Note: due to a lack of plant material, it was not possible to measure leaf starch contents in C. eugenioides .

Similar starch contents for the Coffea arabica allopolyploid and its Coffea canephora diploid parent

Whatever the TRs, there was no significant difference in starch content between C. canephora and C. arabica ( Fig. 2 C). During the day, starch is accumulated in the leaves as an energy source for the coming night. At the time of sampling (ZT = 6 h), the starch contents were very low (2–5% of dry weight) for the two hottest TRs (28–24 and 33–29°C) compared with those obtained (16–22%) for the two coldest TRs (18–14 and 23–19°C). In some plants, cold causes starch accumulation in the chloroplasts. This phenomenon would appear to be very pronounced in the coffee tree.

Level of the phytohormone ABA and its catabolites: original behavior of the allopolyploid under the coldest thermal regime (18–14°C)

The phytohormone ABA and its catabolites played an important role in developmental processes in addition to mediating plant adaptation to stress. The parental diploid species displayed the same pattern, with a strong increase under an intermediate TR (23–19°C; Fig. 2 F). This high ABA content found under this TR was probably the sign of acclimation to an already relatively cold condition for tropical species. Plants suffering from low temperature conditions trigger ABA biosynthesis ( Boudsocq and Laurière 2005 ). The allopolyploid ABA content was comparable with that of the C. canephora parent under the two hottest TRs (33–29 and 28–24°C). Conversely, under the coldest TR (18–14°C), the ABA contents of the two parental species decreased, whereas those of C. arabica reached their highest level.

ABA–glucose ester conjugate (ABAGE) is the most widespread storage form of ABA and considered as a perfect long-distance signal due to the fact that it is translocated to the surrounding tissues. Under the hottest TR (33–29°C), the ABAGE content was similar to those of C. arabica and C. canephora and higher compared with the three other TRs ( Fig. 2 D). For the three hottest TRs (18–14, 23–19 and 28–24°C), ABAGE contents of C. arabica were higher than those of its two parents.

An original behavior of the allopolyploid was found for dihydrophaseic acid (DPA) an ABA oxidation product ( Fig. 2 E). Under the three coldest TRs (18–14, 23–19 and 28–24°C), the DPA content of the two parents was not significantly different and was lower than that of the allopolyploid. As for ABA and ABAGE contents, the allopolyploid DPA content was comparable with that of the C. canephora parent under the hottest TR (33–29°C). It was not possible to know in this study whether DPA accumulation was the result of strong ABA demand, which was catabolized after action, or DPA accumulation for later use in the ABA form.

The allopolyploid plasticity index is lower than that of its parental diploid species

We compared the three Coffea species on the basis of a plasticity index (PI) calculated on the whole transcriptome. By construction, the lower the PI, the less the transcriptome of a species is variable. We found that the PI of the two parental species varied depending on the TR. The C. canephora ’s PI decreases almost linearly from the hottest to the coldest TR. The transcriptome of C. eugenioides is stable at two TRs (28–24 and 23–19°C). Those two TRs mimic the natural conditions for C. eugenioides . The PI increased strongly under the coldest TR (18–14°C), demonstrating a dramatically higher transcriptomic activity. For each TR, the PI of C. arabica was lower than that of its parents ( Fig. 3 ) and appears more stable. These results suggest higher transcriptomic homeostasis of the allopolyploid as compared with the diploid parental species. In addition, we found that the PIs of the two parental species varied depending on the TR and showed a lower index under the two coldest TRs (18–14 and 23–19°C).

Fig. 3

Plasticity index (PI). In order to compare the transcriptome plasticity of each species, a PI of the whole transcriptome for each thermal regime (TR) was calculated.

$$\mathit{P.}\mathit{I}={\displaystyle \sum _{i=1}^n{\left(x_{ij}-\left(y_i\right)\right)}^2/\left({\displaystyle \sum Y_{ji}/a}\right)}$$

where x is the count data (i.e. mean of two replicates) for gene i of species j ; y the mean of counts for gene i for the three species; n = number of genes; Y = the sum of count data for species j ; a = number of species. The lower the PI, the less the transcriptome of a species is variable. Open circles, filled squares and open triangles are for Coffea canephora , Coffea arabica and Coffea eugenioides , respectively. We calculated a standard error by using a bootstrap procedure, removing 20% at a time from each subsample using a circular permutation pattern. For each thermal regime, the plasticity index for C. arabica was lower than that of its parents and appears more stable.

‘Plastic’ and ‘homeostatic’ expression patterns of the three Coffea species

Among the set of expressed genes, fewer ‘homeostatic’ genes compared with ‘plastic’ genes were observed (44% and 56%, respectively; Fig. 4 ).

Fig. 4

Gene expression across three thermal regimes (28–24, 23–19 and 18–14°C). Ten possible expression patterns in Coffea arabica (allotetraploid) relative to its parents ( Coffea canephora and Coffea eugenioides ) across thermal regimes are presented. Each of the 10 patterns involving differential expression is labeled with a Roman numeral and includes a graphical depiction ( Flagel and Wendel 2010 ), where Coffea canephora and Coffea eugenioides are on the outer edges and the polyploid patterns in the middle. Expression values on the same horizontal line indicate statistically equivalent expression, whereas expression values on higher or lower horizontal lines represent statistically significant up- or down-regulation, respectively. For example, the first gene expression pattern (I) represents genes for which no significance was noted in any comparison between the diploids and the allopolyploid whatever the thermal regime considered. Another example is the genes of pattern VIII in which expression in the C. canephora diploid genome was significantly greater than in the C. eugenioides diploid genome, with gene expression in the allopolyploid being significantly different from both but with an intermediate value, whatever the thermal regime considered. Finally, pattern X can represent several patterns (only one graphical depiction is presented here), as it involves genes whose differential expression was variable depending on the thermal regime (i.e. significant differences between thermal regimes and/or significant differences between species under two thermal regimes maximum). The 10 patterns were binned in five categories (color codes), representing ‘no change’ (yellow), ‘ C. eugenioides dominant-like’ (green), ‘ C. canephora dominant-like’ (blue), ‘transgression’ (purple), ‘additivity’ (orange) and ‘not stable’ (gray). The genes of patterns I–IX were considered to be ‘homeostatic’. The genes of pattern X were considered to be ‘plastic’.

Among ‘homeostatic’ genes, the ‘no change’ category accounted for 40% of genes, whereas the categories II–IX (i.e. ‘dominance’, ‘transgressive’ or ‘additivity’) accounted for only 4% of genes. Among the latter, we found 147 genes (VI and VII categories) for which the allopolyploid was in a transgression situation (‘overdominance’).

Among the 56% of genes considered to be ‘plastic’, we determined those genes that had the same expression pattern between species. The hypothesis was that these genes in interaction with the environment were regulated in the same manner. Among them, <1% of the genes had the same expression pattern depending on the TRs ( Supplementary Table S2 ).

‘Non-expressed’ genes present in the three Coffea species whatever the thermal regime

Among studied genes, we focused our attention on genes that we considered as ‘non-expressed’. These were genes whose expression level was close to zero in one species (number of counts <10) and expressed in the two other species (number of counts >250). Only six genes were not expressed in the allopolyploid C. arabica but expressed in the two parental species, notably with three genes encoding Cysteine/Histidine-rich C1 domain family proteins (Cc00_g25320, Cc08_g04790 and Cc08_g04990; Supplementary Table S3 ). These genes are known to be involved in oxidation–reduction processes and are often implicated in defense mechanisms against biotic and abiotic stress. Likewise, only three genes were expressed in the allopolyploid only ( Supplementary Table S3 ).

Studies of key genes in the circadian clock and correlated genes

We investigated some genes linked to the circadian cycle at ZT6, the time at which the greatest divergences in expression were recorded between some diploid species and their polyploid progeny ( Ni et al. 2009 ). These genes are of particular interest, since the circadian clock synchronizes biochemical, metabolic and physiological cycles to environmental changes such as temperatures. It appeared that levels of Late elongated hypocotyl (LHY; Cc02_g39990) and Gigantea (GI; Cc10_g15270) transcripts, two key genes in the circadian clock, were modulated by temperature. For these two genes, we found very strong variations across TRs ( Fig. 5 A, B). The expression of LHY ( Fig. 5 A) was very weak under the two hottest TRs and greatly increased under the highest TR (18–14°C). At this TR, the highest level was reached by the species C. canephora , which was significantly different from the expression level of C. eugenioides and C. arabica . Whatever the TR, the allopolyploid appeared to be in an intermediate situation compared with its two parents. The increase in the level of LHY mRNA with a temperature decrease was consistent with the observations of Gould et al. (2006) . Chlorophyllide a oxygenase (CAO; Cc10_g11980), which is the key enzyme in Chl b biosynthesis ( Mochizuki et al. 2010 ), is known to respond strongly to circadian signals. Indeed, we found a strong correlation with LHY ( r > 0.94; Supplementary Table S4 ). The antenna size, an important factor for the efficiency of photosynthesis, is closely associated with the biosynthesis of Chl b . The increase in CAO expression in C. canephora under the coldest TR seemed to correspond to an increase in Chl degradation to reduce the phototoxic effects of Chls released from the photosystems ( Figs. 1 B, 5D ).

Fig. 5

Expression variations for the three Coffea species of LHY, GI, CAO, AAO and catalase 3 genes and the mRNA expression level ratio for GI/LHY depending on the thermal regime. (A) Late elongated hypocotyl (LHY; Cc02_g39990), (B) Gigantea (GI; Cc10_g15270), (C) GI/LHY ratio, (D) chlorophyllide a oxygenase (CAO; Cc10_g11980), (E) ascorbate oxidase (AAO; Cc07_g15610) and (F) catalase 3 (Cc10_g00570). Open circles, filled squares and open triangles are for Coffea canephora , Coffea arabica and Coffea eugenioides , respectively. Asterisks indicate significant differences (* P < 0.05; ** P < 0.01) between the particular parent and the allopolyploid for each thermal regime.

As regards the GI gene, we found strong down-regulated expression for the three species under the intermediate TR 28–24°C ( Fig. 5 B). Compared with its parents, C. arabica displayed identical levels or was in an additive situation. Under the coldest TRs (18–14 and 23–19°C), the allopolyploid displayed an expression level similar to that of C. canephora . Under the hottest TR (33–29°C), the expression level of C. canephora was significantly up-regulated compared with that of the allopolyploid. The pattern of GI expression was highly significantly and linearly correlated ( Supplementary Table S5 ) to the patterns of Adagio (a Flavin-binding Kelch repeat, F BOX protein; Cc11_g07690), a member of the Zeitlupe (ZTL) protein family.

A dynamic balance between LHY and GI is necessary for effective temperature compensation at high temperatures ( Gould et al. 2006 ). These authors have shown that the levels of GI could be modulated by temperature and that the GI gene plays a critical role in increasing the temperature range permissive for rhythmicity ( Gould et al. 2006 ). The coupling between GI and LHY can be visualized in the form of the ratio of GI transcripts to LHY transcripts ( Fig. 5 C). Under the highest TR (33–29°C), the very low LHY expression levels appeared to be counterbalanced by an increase in GI, as observed by Gould et al. (2006) and Ni et al. (2009) for diploid species as well as for allopolyploids. We found that, under the highest TR (33–29°C), the ratios were very high for C. canephora and the allopolyploid. There appeared to be a strong correlation ( r > 0.92) between the level of ABAGE and the GI/LHY ratio. Some strong correlations were seen between the GI/LHY ratio and genes which belonged to the pentatricopeptide repeat (PPR) gene family (Cc00_g20160, Cc01_g13270, Cc01_g13440, Cc03_g04950, Cc05_g05570, Cc06_g06070, Cc06_g08820, Cc07_g16390 and Cc08_g10650; Supplementary Table S6 ).

Genes involved in photosynthesis, respiration and reactive oxygen species

We examined the response of allopolyploid genes involved in the biosynthesis of photosynthetic pigments and in the ‘light reactions’ of photosynthesis ( Supplementary Table S7 ), respiration ( Supplementary Table S8 ) and ROS detoxification/production ( Supplementary Table S9 ), since the performance of these phenomena is dependent upon ambient temperatures. No gene was found for which C. arabica appeared to be in a situation of overdominance compared with its two parents under all the TRs. A few genes for which we found high expression differentials between the parental species and the allopolyploid were studied further.

Under the hottest TR (33–29°C), there was strong down-regulated expression of C. arabica compared with its parents for the genes of the light-harvesting complex (LHCII; Cc04_g16410) and Chl a-b binding protein (CAB; Cc10_g00140, Cc05_g12720, Cc09_g09020, Cc05_g09650 and Cc09_g09030) associated with PSII (Cc06_g11950, Cc10_g11890 and Cc10_g11890; Supplementary Table S7 ).

There were few genes involved in respiration whose level of expression was strongly affected by allopolyploidy. It was particularly under the coldest TRs (18–14 and 23–19°C) that C. arabica stood out from its C. canephora parent ( Supplementary Table S8 ). For example, up-regulated expression of the Calvin cycle gene sedoheptulose-1,7-bisphosphatase (Cc02_g06960) was observed under these two TRs for the C. arabica / C. canephora comparison.

For all the genes linked to ROS ( Supplementary Table S9 ), there were several genes that showed down- and/or up-regulated expression in the allopolyploid compared with its parents. Many of them are known for their role as ROS scavengers (catalase, Cc10_g00570; thioredoxin, Cc06_g21570, Cc10_g02320, Cc08_g11260, Cc01_g14690 and Cc09_g03800; peroxidase, Cc02_g30380, Cc00_g17550, Cc07_g02500, Cc07_g06870 and Cc05_g08480; and glutathione, Cc08_g14620). The allopolyploid showed 17-fold up-regulated expression compared with C. canephora and 11-fold compared with C. eugenioides . Conversely, the allopolyploid showed 27-fold down-regulated expression compared with C. canephora and 29-fold compared with C. eugenioides . Under the TRs considered to be the most stressful (18–14 and 33–29°C) for the three species studied, a corresponding up-regulation of genes involved in the detoxification of ROS was expected. Indeed, abiotic stress results in an imbalance of cell redox status due to the overproduction of oxidative radicals and can cause damage to a wide range of cellular components, such as the photosynthetic apparatus and various other components, thereby hindering metabolic activities and affecting plant growth by limiting metabolic flux activities ( Sairam and Tyagi 2004 ). Ascorbate oxidase ( AAO , Cc07_g15610) encodes an apoplastic enzyme that catalyzes the reversible oxidation of ascorbate to 2-dehydro-ascorbate with the concomitant reduction of molecular oxygen to water. Here we observed that the expression level of C. arabica was the same as that of C. canephora under the hottest TR (33–29°C; Fig. 5 E) and comparable with that of its two parents under intermediate TRs (28–24 and 23–19°C). On the other hand, under the coldest TR (14–18°C), C. arabica had a significantly higher expression level than that of its parents.

For catalase 3, which displayed high differential expression (Cc10_g00570; Fig. 5 F), the allopolyploid appeared to be in an additive situation whatever the TR compared with its two parents. In addition, the expression of this gene was stable depending on the TRs.

Finally, polyphenol oxidases (PPOs) catalyze the oxidation of various phenolic substrates to quinones. PPOs are implicated in defense against insects, plant pathogens and mechanical damage. In coffee, the constitutive levels of PPO activity are relatively higher than in other species ( Mazzafera and Robinson 2000 ). Very high overexpression of three PPO genes was seen in C. canephora under the coldest TR (14–18°C; PPO: Cc05_g10250, Cc00_g35890 and Cc05_g10310; Fig. 6 A–C). Three ascorbate peroxidases (APX family) were studied ( Supplementary Table S9 ). Only APX2 (Cc01_g11800), which is targeted to the cytosol, displayed a high expression level in negative overdominance of the allopolyploid compared with its C. eugenioides parent and to a lesser degree compared with C. canephora ( Fig. 6 D). The other two genes did not display any differences between species.

Fig. 6

Expression variations for the three Coffea species for PPO, APX2 and FAD7/8 genes depending on the thermal regime. (A) Polyphenol oxidase chloroplastic (PPO; Cc05_g10250), (B) polyphenol oxidase chloroplastic (PPO; Cc00_g35890), (C) polyphenol oxidase chloroplastic (PPO; Cc05_g10310), (D) l -ascorbate peroxidase 2 cytosolic (APX2; Cc01_g11800) and (E) fatty acid desaturases (FAD7/8; Cc02_g06400). Open circles, filled squares and open triangles are for Coffea canephora , Coffea arabica and Coffea eugenioides , respectively. Asterisks indicate significant differences (* P < 0.05; ** P < 0.01) between the particular parent and the allopolyploid for each thermal regime.

Gene expression patterns highly correlated with metabolites

The contents of metabolites for which fluctuations are observed between species have been correlated with gene expression. For linolenic acid (C18:3), which results from the desaturation of linoleic acid (C18:2) by omega-3 fatty acid desaturases (FADs), two enzyme orthologs of both AtFAD3 and AtFAD7/8 contribute to C18:2 desaturation and were identified in coffee ( Denoeud et al. 2014 ). The transcript CcFAD3 (Cc01g14330) was very poorly expressed in coffee leaves (data not shown), whereas CcFAD7/8 (Cc02g06400) exhibited high transcript levels in the three species, suggesting a major role for this gene in the leaf C18:3 content. The CcFAD7/8 expression level displayed differential expression between species ( Fig. 6 E) We observed that the allopolyploid transcript level was at the same level as that of C. canephora under the hottest TRs (28–24 and 33–29°C) and significantly higher than that of C. eugenioides , but at a higher level than that of C. canephora under the coldest TRs (18–14 and 23–19°C). For other metabolites such as sucrose, stacchyose, rafinose, ABA and its catabolites, we did not observe any significant correlation with upstream genes in biosynthetic pathways. However, these metabolites are strongly coupled to the expression of genes that are part of signaling and regulatory systems ( Supplementary Table S10 ). For instance, MYB44 (Cc07_g06450) is highly and negatively correlated ( r = –0.864) with the concentrations in ABAGE. MYB transcription factors are key components in ABA signaling ( Yamaguchi-Shinozaki et al. 2006 ). Jung et al. (2008) showed that ABA induced MYB44 transcript accumulation within 30 min. The gene was also activated under various abiotic stresses, and overexpression of MYB44 enhances stomatal closure to confer abiotic stress tolerance in transgenic Arabidopsis. The ABAGE contents are also highly correlated with ERF4 (Cc05_g07590), known as a transcriptional repressor capable of modulating ethylene and ABA responses ( Yang et al. 2005 ). Increased sucrose levels lead to increased levels of RFOs and anthocyanins ( Teng et al. 2005 ). Here the high correlation between the raffinose contents and the expression level of the malonyl-coenzymeA:anthocyanin 3- O -glucoside-6′′- O -malonyltransferase (Cc05_g05980; Supplementary Table S10 ) seems to confirm the synergistic interaction of sugars and phenolic compounds forming part of an integrated redox system that contributes to stress tolerance.

Raffinose is negatively correlated with calmodulin-CIP111 (Cc06_g20150; Supplementary Table S10 ). The calmodulin family is a major class of calcium sensor proteins through the regulation of numerous target proteins. The induction of calmodulin genes by abiotic stimuli has been described, and genetic studies indicate positive as well as negative effects of calmodulin-regulated pathways on stress responses ( Ranty et al. 2006 ).

Discussion

The purpose of our study was to measure the performance of a recently formed allopolyploid species and explain, on a molecular basis, the differences in performance observed between the allopolyploid and its parental diploid species. In order to observe the differences between species, plants were subjected to a range of thermal regimes (TRs) for long enough to enable important changes to be measured, involving changes at the phenotype, transcriptome and physiological levels.

Growth performance: the allopolyploid Coffea arabica was more homeostatic

As already detailed in the Introduction, C. canephora tolerates rather high growth temperatures: 22–26°C on average ( DaMatta and Ramalho 2006 ). In our study, C. canephora displayed very good adaptation to the hottest TRs, with an optimum occurring under the 28–24°C TR. These growth performances tallied with that observed under natural conditions. Under the coldest TR (18–14°C), the C. canephora growth rate collapsed and photobleaching symptoms were seen on most of the seedlings. They showed real physiological disorders probably reflecting some extensive damage to the chloroplasts. These observations corresponded to what is commonly observed for cultivated C. canephora varieties ( Partelli et al. 2009 ). For its part, C. eugenioides grows naturally at rather high tropical elevations, which consequently display a rather cool climate (18–23°C according to Davis et al. 2006 ). In our study, C. eugenioides displayed very good adaptation to the coldest TRs, with an optimum occurring under the 23–19°C TR. When the TR was very hot (33–29°C), there was a very low growth rate associated with high mortality, reflecting poor adaptation of the species. Our study therefore confirmed that C. eugenioides is a species that is rather well adapted to cold high altitude climates. As regards the allopolyploid species C. arabica , its behavior was similar to that of the parental species best adapted to the TR considered. Under the hottest TR (33–29°C), its growth rate was identical to that of C. canephora , and under the coldest TR (18–14°C) identical to that of C. eugenioides . Unlike its parental species, C. arabica did not display any major variations at the extreme TR levels, but a stability along the thermal range. This stability reflected phenotypic homeostasis, while the variations seen in the parents reflected greater phenotypic plasticity. Our observations for C. arabica corresponded to what has been seen under growing conditions, as this species is currently grown at average temperatures between 18 and 26°C. It has already been found that C. canephora shows significant decreases in growth under the 18–14°C TR when compared with C. arabica ( Ramalho et al. 2003 ). On the other hand, at temperatures >35°C, the coffee plants of both species seem to possess the same tolerance under growing conditions ( DaMatta et al. 2008 ). In fact, under intensive management conditions it is possible to grow Arabica coffee plantations in marginal regions with average temperature as high as 24–25°C ( DaMatta et al. 2008 ).

The metabolic levels of allopolyploid Coffea arabica are remarkably similar to those of Coffea canephora under the hottest thermal regime (33–29°C)

Under the highest TR (33–29°C), the behavior of the allopolyploid was remarkably comparable with that of its C. canephora parent. This phenomenon was found for all the metabolite contents measured in this study. As was observed by DaMatta et al. ( 2001 , 2006 ), the temperature needed to obtain A_max (photosynthetic capacity, obtained under light- and CO ₂ -saturating conditions and optimal temperature) was as high as 35°C for both Arabica and Robusta coffee. Under the intermediate 28–24°C TR, C. arabica often displayed a behavior similar to that of C. canephora , except for the total sugars, sucrose, stachyose and DPA contents, where its behavior seemed to be original, and for raffinose, where its performance was similar to that of the C. eugenioides parent. Under the two coldest TRs (18–14 and 23–19°C), the behavior of C. arabica differed depending on the measurement taken. All the same, it was noted for the coldest TR (18–14°C) that the behavior of C. arabica was often original compared with the two parental species, as for the sugars, or similar to that of C. eugenioides , as was the case for the Chl and linolenic acid contents. As was explained in Damatta et al. (2006) , the effect of low positive temperature was studied on C. arabica and C. canephora cultivars ( Quartin et al. 2000 , Campos et al. 2003 , Ramalho et al. 2003 ), and it was shown that C. arabica cultivars exhibited better cold acclimation ability than the other species studied. On the other hand, Campos et al. (2003) showed that changes in lipids would largely explain the better cold tolerance of C. arabica cultivars since such tolerance has been related to a much lower fatty acid saturation level (mainly C16:0 and C18:0) in phosphatidylglycerol than in cold-sensitive cultivars ( Damatta et al. 2006 ).

The allopolyploid’s overall transcriptional response to the thermal regime was more homeostatic

We showed by calculating a PI onto the whole transcriptome that the allopolyploid was more homoeostatic for global gene expression than its parental diploids species. This means that the gene expression levels of the allopolyploid were less variable in relation to the environment than those of the diploid parents. On the other hand, the latter displayed greater plasticity. The possession of two complete subgenomes from different parental species seems to promote homeostasis in allopolyploids compared with diploid species. A similar result was recently published by Coate et al. (2013) , in an allopolyploid Glycine dolichocarpa. However, it remains to be determined whether expression homeostasis for a large number of genes may provide a molecular basis to explain adaptation to variable ambient temperature ranges or, conversely, whether it is the consequence of prior acclimation.

The majority of gene expression in strong interaction with the thermal regime

To investigate the effects of genomic merger and doubling, transcriptomic divergence between parents and synthetic nascent allopolyploids, or natural allopolyploids, has been assessed in several recent studies using genome-wide approaches to measure the deviation from additivity ( Hegarty et al. 2006 , Wang et al. 2006 , Hegarty et al. 2008 , Chaudhary et al. 2009 ). In our study, we applied the same partitioning rule as Rapp et al. (2009) , adding on this additional condition: in order for a gene to be classed as ‘additive’, ‘dominant’ or ‘overdominant’, it had to be so for all the TRs. Around 40% of genes were considered as ‘homeostatic’, because they displayed no significance (‘no change’ category) in any comparison between the diploids and the allopolyploid whatever the TR considered. This large proportion showed to what degree the three species shared similar expression levels under TRs that were nonetheless highly contrasting for the species studied. Only 4% of genes were differentially expressed between species and for which the allopolyploid was systematically positioned in a situation of dominance, overdominance or additivity. Among the overdominance genes, we noted for future in-depth study the presence of the Hasty 1 gene, a member of the importin-β family in Arabidopsis and known to regulate the nucleocytoplasmic transport of molecules as micro RNAs and, for this reason, involved in several different morphogenetic pathways ( Park et al. 2005 ). For the remaining 56% of expressed genes, ‘plastic’ expression patterns were found, i.e. in strong interaction with the TR. These so-called ‘plastic’ genes had variable relations of dominance, overdominance and additivity between TRs. This observation backed up our earlier results whereby we showed that some genes in an additive situation for one TR became, for example, ‘ C. canephora -like dominance’ for another TR and more globally that the transcriptome divergence between C. canephora and C. arabica was modulated by the TR ( Bardil et al. 2011 ).

Until recently, we believed that there existed a close link between expression level dominance and homeolog expression bias in coffee as has been shown in cotton ( Yoo et al. 2013 ). Indeed, it was considered that ‘transcriptomic shock’ ( Hegarty et al. 2006 , Buggs et al. 2011 ) generated by the merging of two different genomes resulted in massive novel patterns of homeolog activation and repression, through massive transactivation and repression due to the divergence in parental regulatory mechanisms that become reunited in allopolyploids ( Yoo et al. 2013 ). The results obtained by our team with C. arabica ( Combes et al. 2012 , Combes et al. 2013 ) showed that the subgenome contributions to the transcriptome appeared close to balanced and marginally altered by the growing temperatures. We have concluded that C. arabica ’s ability to tolerate a broader range of growing temperatures did not result from massively differential use of homeologs. In the present study we showed a very limited number of ‘non-expressed genes’. The large proportion of ‘plastic’ genes on the one hand and the low homeolog expression bias on the other hand suggests that differential gene expression is probably mainly dependent on a small number of regulators (transcription factors, epigenetic mechanisms, transposons and post-transcriptional regulation), which are themselves probably highly dependent on the ambient temperatures. Among the genes expressed only in the allopolyploid, we found, for example, the Cc00_g14760 gene, which is a mitogen-activated protein kinase (MAPK) that resembles AT1G51660. The MAPKs play important roles at the intersection between environmental stress responses and internal growth programs. They are involved in plant immunity, cell death regulation or development of stomata or epidermal spores ( Takemoto et al. 2005 , Kim et al. 2012 ). A second gene of interest for future studies was Cc03_g12020, very similar to AT3G42170, which is a transposase, called DAYSLEEPER, that is essential for normal plant growth ( Bundock and Hooykaas 2005 ). This transposase-like gene can regulate global gene expression.

No massive gene reprogramming is linked to photosynthesis and respiration of the allopolyploid Coffea arabica

We showed that the allopolyploid state did not lead to massive reprogramming of allopolyploid genes involved in photosynthesis light reaction, carbon fixation and photorespiration, because the gene expression differences observed were few in number and were not generally very large. For three genes encoding ferredoxin (Fd; Cc02_g28520, Cc06_g19130 and Cc07_g10820), we found up-regulated expression of the allopolyploid compared with C. canephora under the coldest TR (18–14°C). In chloroplasts, Fd accepts electrons from PSI and transfers them to the enzyme Fd:NADP ⁺ oxidoreductase for the photoreduction of NADP ⁺ ( Ramirez et al. 2013 ). It has been demonstrated that Fd activity is one of the rate-limiting steps in photosynthesis, and controls the balance between the demand for redox equivalents and photosynthetic activity under a wide range of environmental conditions ( Hajirezaei et al. 2002 ). Fd also participates in fatty acid synthesis.

The Calvin cycle sedoheptulose-1,7-bisphosphatase gene expression level did not differ between the allopolyploid and C. eugenioides , but was significantly higher than in C. canephora under the coldest TR (18–14°C). In transgenic plants of Arabidopsis, overexpression of this gene is reflected in an increase in biomass and higher sucrose and cumulated starch levels during the photoperiod ( Lefebvre et al . 2005 ). We found a higher sucrose content under the coldest TR (18–14°C) in the allopolyploid.

However, for the majority of genes, C. arabica usually had an expression level that was ‘additive’ or identical (‘no change’ category) compared with its two parents. In general, there were few notable differences between the species, and the reasonable hypothesis could be put forward that the rate of carbon assimilation differed little from one species to the other under the intermediate TRs (23–19 and 28–24°C), in contrast to observations made in other species ( Ni et al. 2009 , Coate et al. 2013 ).

No major alteration of the circadian clock genes occurs between the allopolyploid Coffea arabica and its parental diploid species

Ni et al. (2009) showed in some allopolyploid neosynthetics of Arabidopsis that an alteration of the circadian rhythm was directly linked to increased expression of downstream genes and metabolic pathways, leading to increased amounts of Chls, starch and sugars during vegetative growth. Here, we only studied the circadian clock at ZT6, which was no doubt not enough to confirm that the circadian clock matched the environment in the same way in the three species. For two major genes of the circadian cycle (LHY and GI), the allopolyploid displays significant differences from one or both parents only under the coldest and hottest TRs (18–14 and 33–29°C).

Among the genes linearly correlated with LHY expression, there was also a serine/arginine-rich splicing factor. These proteins are emerging as one of the master regulators of gene expression ( Long and Caceres 2009 ). Recent data suggest that alternative splicing may be involved in the regulation of clock temperature compensation ( Filichkin and Mockler 2012 ). Expression of the Cc02t38220.1 gene, a ubiquitin-conjugating enzyme E2s, was itself also highly correlated to that of the LHY gene ( r > 0.97). Many intracellular proteins become covalently modified with ubiquitin or ubiquitin-like proteins, and these modifications determine protein turnover, regulation and molecular function, providing a complex regulation of the proteome ( Wijk and Timmers 2010 ). The expression of LHY is also strongly regulated with that of histone deacetylase 19 (HDA 19). This enzyme catalyzes deacetylation at several loci in response to abiotic stress ( Bilichak and Kovalchuk 2013 ).

Among the genes correlated to GI, we noted several representants of the ZTL family. Baudry et al. (2010) indicated that all members of the ZTL protein family are engaged in proline-rich protein 5 (PRR5) and Timing of cab expression 1 (TOC1) degradation, two essential components of the oscillator. Also worth noting were the strong correlations between GI and the CABs (CAB21, CAB36 and CAB13). CAB genes encode proteins which modify the plane of orientation of chlorophylls. Nine genes belonging to the PPR family were highly correlated with the GI/LHY ratio. PPR proteins ( Lurin et al. 2004 ) have a range of essential functions in post-transcriptional processes (including RNA editing, RNA splicing, RNA cleavage and translation) within mitochondria and chloroplasts. It has been shown that some PPRs regulate ROS homeostasis in mitochondria by interacting with genes such as alternative oxidase ( Laluk et al. 2011 ).

We concluded that the alteration of the circadian clock presented as a rule in other allopolyploid species seems to depend mainly on the ambient cues for C. arabica . Heterosis would therefore appear not to be the compulsory prerogative of allopolyploid species, and allopolyploids do not necessarily display a ‘fixed’ alteration of the circadian cycle, in contrast to what was assumed by Chen (2010) . Coffea arabica is a counter-example, which is doubtless not alone among allopolyploid species.

Greater capacity of the allopolyploid Coffea arabica for membrane fluidity and reactive oxygen species homeostasis

One of the main causes we chose to explain the greater homeostasis of allopolyploid growth across the contrasting TRs was the capacity of the allopolyploid to maintain better membrane fluidity and integrity. In fact, an important result of our study concerned fatty acids, ROS and the RFOs. Long-term exposure to the range of TRs in our experiment promoted changes in the balance of the fatty acid composition of the three species considered. In that respect, the original behavior of the allopolyploid needs to be mentioned. Indeed, under the two hottest TRs (28–24 and 33–29°C), the linolenic acid content of the allopolyploid leaves was similar to that of C. canephora , whereas under the two coldest TRs (18–14 and 23–19°C), their composition was similar to that of C. eugenioides . Hence, we are seeing great plasticity of the allopolyploid depending on the thermal range. Increases in unsaturated fatty acids in cellular membranes, especially at the beginning of acclimation, have been reported in coffee plants ( Campos et al. 2003 ), and recently Scoti-Campos et al. (2013) showed that C. arabica displayed an earlier acclimation response than C. canephora at low temperatures. Under the highest TR (33–29°C), it is likely that the membrane lipid composition of C. eugenioides was a highly unfavorable factor for good acclimation since the majority of the plants died. In addition, we found that stacchyose and raffinose contents were higher in the allopolyploid than in its two parents for the coldest TRs. RFOs act as stabilizers of cellular membranes but also as ROS scavengers ( Nishizawa et al. 2008 ). Dos Santos et al. (2011) previously reported an accumulation of RFOs in coffee leaves in response to different types of stress (heat and water deficit). It is also known that exogenous supplies of ABA lead to an increase in stachyose ( De Block et al. 2005 ). Indeed, we found higher contents or ABA and its catabolites in the allopolyploid under the coldest TR that could explain the accumulation of RFOs.

The coldest TR (18–14°C) was definitely stressful for C. canephora . The poor growth and substantial loss of Chl in that species would seem to be due to its inability to prevent and manage oxidative stress. In coffee, the best known action process is based on its capacity to oxidize phenolic compounds when tissue is damaged ( Melo et al. 2006 ). In tomato, PPO overexpression increased photo-oxidative damage during water stress ( Thipyapong et al. 2004 ). The high overexpression of PPOs observed in this study is an indicator of the degree of stress to which this species was subjected under the coldest TR. PPOs are powerful ROS producers and doubtless reflected a disorganization of membranes in C. canephora under the coldest TR. However, it was not possible to conclude whether the strong overexpression of these three genes in C. canephora was a consequence of stress or an inducer of stress.

As we were unable to study the transcriptome of C. eugenioides under the hottest TR (33–29°C), we can merely put forward hypotheses. It is likely that the very high mortality and the very low growth rates for C. eugenioides under this TR were due to some detrimental effects of heat on the photosynthetic apparatus associated with a poorly adapted membrane lipid composition and to the inability of this species to re-establish cellular redox homeostasis.

In C. arabica , redox activity appeared more stable than that of the diploid parents. We put forward the hypothesis that the increase in AAO expression found in the allopolyploid under the coldest TRs reflected better redox homeostasis than in its two parents. AAO appears highly expressed in C. arabica and could be a first line of defense against the oxidative stress. AAO, which is a key regulator of extracellular redox status ( Noctor 2006 ), also has a role in the perception of the environment or stress responses ( Pignocchi et al. 2003 ). Coate et al. (2013) showed that under chronic excess light the capacity for non-photochemical quenching was higher in a recently formed natural allopolyploid ( Glycine dolichocarpa ) than in its diploid parents ( G. tomentella and G. syndetika ). Worth noting is that expression of the WRKY33 gene was highly correlated to that of AAO ( r > 0.91). Li et al. (2011) showed that this transcription factor is involved in plant thermotolerance. Regarding the APX genes, it is known that the multiple isoforms throughout the chloroplast, mitochondrion, peroxisome and cytosol are differentially regulated in a tissue-dependent manner during oxidative stress responses ( Rossel et al. 2007 ), so it is not surprising that only APX2 was differentially expressed. The increase in APX2 under an intermediate TR (23–19°C) may have corresponded to the increase in ABA observed since the APX2 promoter contains ABA-responsive elements ( Rossel et al. 2007 ).

Transcriptional homeostasis of the allopolyploid Coffea arabica : better capacity to sense temperatures

Here, we put forward the hypothesis that the primary temperature sensing ( Mittler et al. 2012 ) of the allopolyploid is more sensitive (i.e. more efficient) than that of the parental species. That greater sensitivity would seem to trigger quicker reprogramming of the transcriptome, proteome and metabolome, in a cascade, in order to adapt to contrasting TRs and thermal stress. Strong correlations between hormones, sugars, linolenic acid and transcription factors suggest that in turn variations of these metabolites play a crucial role in cell signaling cascades. On the other hand, in the diploid parental species, the signals could seem to occur later and trigger more massive, but less effective programming. This could explain the lower homeostasis of their transcriptome. As plasma membranes contain combinations of glycosphingolipids and protein receptors, some optimum modifications in membrane fluidity in the allopolyploid are apt to explain how these sensors might be affected more rapidly.

Several questions therefore arise for future research. The greater homeostasis in the allopolyploid might lie in its ability to organize changes in membrane fluidity. Testing this hypothesis will call for an experimental design making it possible to study expression kinetics.

Conclusion

To date, many allopolyploids have seemed remarkable because they introduced a very great and obvious break from their parents, pointing towards greater phenotypic plasticity as an important determinant in species invasiveness ( Davidson et al. 2011 ). However greater homeostasis along environmental gradients has also been suggested to explain the evolutive success of polyploids ( Bretagnolle and Thompson 2001 ).

Here, we conclude that better homeostasis of the allopolyploid C. arabica transcriptome affords greater phenotypic homeostasis when faced with environments that are unsuited to the parental species and leads to better adaptive ability in the allopolyploid species. Our results lead to a new vision of the potential of the C. arabica allopolyploid compared with its parental diploid species. In particular, they help to explain how a species from a tiny Sudano-Ethiopian ecological niche has been grown since recent times (around two centuries) throughout the intertropical belt in different climates, and all this from a weak genetic diversity.

In view of predicted global warming, understanding how the membrane lipid saturation and protection regulation system is modified in the allopolyploid and why that affords it an intrinsic adaptation advantage is a promising area to develop homeostatic varieties that are more capable of coping with adverse environmental conditions.

Materials and Methods

Plant material and cultivation in phytotron chambers

Fresh mature seeds of C. arabica (L.) and C. canephora (Pierre) were provided by the ‘Nicafrance’ Experimental Station in La Cumplida (Matagalpa, Nicaragua). Coffea arabica was represented by the cultivar ‘Caturra’. Coffea canephora was represented by the cultivar ‘Nemaya’, derived from a cross between two wild Congolese genotypes ( Bertrand et al. 2000 ). Finally, mature C. eugenioides seeds were provided by the ‘Coffee Research Foundation’, from cultivated trees in a collection and originating from the Mount Elgon forest in Kenya (1°8′N 34°33′E).

The coffee seedlings were grown in a glasshouse under natural daylight, at a constant temperature of 24°C. After 5 weeks, the dry weights (103°C for 24 h) were determined for 45 seedlings of each species. This dry weight was used as a reference to account later for the growth of seedlings subjected to different temperature conditions. For each species, 120 seedlings were shared out into 12 mini-glasshouses. With a view to uniformity, the three species studied were represented in each of the mini-glasshouses by 10 seedlings per species, i.e. 30 seedlings per mini-glasshouse. Four day/night temperature conditions (called TRs) were tested: 18–14, 23–19, 28–24 and 33–29°C.

For each TR, three mini-glasshouses were transferred per phytotron (CRYONEXT, model RTH 1200L). The following parameters were fixed for each phytotron: photoperiod (12 h/12 h light/dark), humidity (80%) and luminosity (600 µmol m ⁻² s ⁻¹ ). After 11 weeks in phytotrons, five seedlings were sampled per species and per mini-greenhouse. The dry weights (103°C for 24 h) were determined to establish growth ratios (expressed as a pourcentage of dry weight) compared with the reference. On the remaining seedlings, leaf samples were taken for RNA extraction and chemical analyses. For RNA extraction, leaves of the terminal young pairs were collected at subjective midday, corresponding to ZT6, and then flash-frozen in nitrogen and stored at 80°C until RNA extraction. At the same time, for the chemical analyses, leaves of the terminal young pairs were sampled, flash-frozen in nitrogen and then freeze-dried. During sampling operations, three biological replicates were performed for chemical analysis and two for RNAseq.

Chemical analyses

All analyses were carried out in triplicate, from three different extractions, using a completely random experimental design.

For determination of Chl content, we used the Mc Kinney (1941) and Holden (1976) methods. A standard range was established using a pure extract of Chl a and b from spinach (Sigma-Aldrich C5753-1MG and C5878-1MG). The Chl content was calculated using spectrophotometric absorbance at light wavelengths of 603, 645 and 663 nm.

Total lipids were extracted from 100 mg samples of freeze-dried powder using a modified Folch method as described previously by Laffargue et al. (2007) . Fatty acid methyl esters (FAMEs) were prepared according to the ISO-5509 standard, and gas chromatography analyses of FAMEs were performed as previously described by Laffargue et al. (2007) .

Sugars were extracted from 20 mg samples of freeze-dried powder and measured by high-performance anion exchange chromatography coupled with pulsed amperometric detection (Dionex Chromatography Co.) as described by Dussert et al. (2006) .

Starch contents were determined on 800–1,000 mg of dry powder, using a starch enzymatic method (NF V18-121), by InVivo Labs (Vannes, France). After elimination of soluble sugars and the soluble products of starch degradation, the residue was dispersed using sodium hydroxide. The starch was then hydrolyzed into glucose units with amyloglucosidase (D-Glucose-HK enzyme kit, Megazyme). Lastly, the glucose obtained was quantified by spectrophotometry using a hexokinase system. Due to a lack of plant material, it was not possible to measure leaf starch contents in C. eugenioides .

ABA and its catabolite profiling analysis were done by the National Research Council Canada. ABA and its metabolites were quantified by ultra-performance liquid chromatography–elecrospray ionization–tandem mass spectrometry as described in detail by Chiwocha et al. (2003) .

RNA extraction

Total RNA was extracted with a protocol modified from Morcillo et al. (2006) based on Corre et al. (1996) . For each sample, approximately 1 g of leaves was used, except for the RNA extraction of C. eugenioides leaves under the hottest TR (33–29°C; not enough biological samples). The assay of total RNA was carried out by measuring the optical density of a solution of 1 µl at 260 nm using a spectrophotometer adapted to small volumes (Nanodrop, ND, Labtech).

Validation of RNAseq by qRT-PCR analysis

In RNAseq studies, RT-PCR is usually used for validation.

Based on bibliographic data, we targeted five key genes of the circadian cycle [Circadian clock-associated 1 (CCA1), TOC1, GI, ZTL and PRR5) and four housekeeping genes [ubiquitin (UBQ), 60S ribosomal protein L7 (RP17), 14-3-3 protein and clathrin-associated protein (AP47)] as an internal control. The expression of those genes was measured in the four thermal regimes and for the three species. The protocol is given in Supplementary Table S11 .

Normalization was performed as previously described by Salmona et al. (2008) and Joët et al. (2009) . Briefly, normalization was processed by using the mean expression (C _t value) for all genes under study. Relative expression rates of target gene transcripts were then calculated using the following formula:

$${\text{E}}^{(-}{}^{\Delta \text{Ct})}$$

where Δ _Ct = C _t target gene − C _t mean of four housekeeping genes and E = efficiency value for each gene and for each sample.

As in Salmona et al. (2008) and Joët et al. (2009) , the expression stability of four putative housekeeping genes selected from the literature on Coffea species ( Barsalobres-Cavallari et al. 2009 , Cruz et al. 2009 ) was determined to obtain an accurate internal normalizer.

We showed that the fold changes estimated from RNAseq had high correlation with that from the qRT-PCR ( r > 0.80; Supplementary Table S11 ).

RNA sequencing and bioinformatics

RNAseq was carried out by the MGX platform (Montpellier GenomiX, www.mgx.cnrs.fr/ ). RNAseq libraries were constructed with the TruSeq RNA sample preparation (Low throughput protocol) kit from Illumina. A 1 µg aliquot of total RNA was used for the construction of the libraries. The first step in the workflow involves purifying the poly(A)-containing mRNA molecules using poly(T) oligo-attached magnetic beads. Following purification, the mRNA is fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments are copied into first-strand cDNA using SuperScript II reverse transcriptase and random hexamer primers. This was followed by second-strand cDNA synthesis. These cDNA fragments then go through an end repair process, the addition of a single ‘A’ base and subsequent ligation of the adaptor. The products are then purified and enriched with 15 cycles of PCR. The final cDNA libraries were validated with a DNA 1000 Labchip on a Bioanalyzer (Agilent) and quantified with a KAPA qPCR kit. For each sequencing lane, four or six libraries were pooled in equal proportions, denatured with NaOH and diluted to 6.5 pM before clustering. Clustering and 100 nucleotide single read sequencing were performed according to the manufacturer’s instructions.

Image analyses and basecalling were performed using the HiSeq Control Software and Real-Time Analysis component (Illumina). The quality of the data was assessed using fastqc from the Babraham Institute ( http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ ) and the Illumina software SAV (Sequence Analysis Viewer). Using the Illumina-supplied algorithm CASAVA (1.8.2 version), which is designed to be particularly efficient for 74 bp reads, 72 bp reads were generated from libraries derived from RNA samples. The 72 bp reads were aligned using BWA (Burrows–Wheeler Aligner; Li and Durbin 2009 ) against the C. canephora 25,574 protein-coding gene models used as the reference transcriptome ( Deneoud et al. 2014 ). A maximum of two mismatched nucleotides (including gaps) was allowed between the read and the reference. Unmapped and multiposition matched reads were excluded from analyses, whereas uniquely mapped reads were counted for each gene model using samtools for subsequent differential expression analysis by DESeq ( Anders and Huber 2010 ). By these criteria, 66–69% of reads mapped uniquely to a genomic location ( Supplementary Table S12 ).

Differential expression analysis

The RPKM (reads per kilobase per million mapped reads) data are provided in Supplementaty Table S12 . However, to compare the expression levels between the species and between the growth TRs, we used the DESeq methodology ( Anders and Huber 2010 ) which is more accurate than RPKM data ( Dillies et al. 2013 ; Supplementary Table S13 ). The DESeq package is particularly adapted to experiments with biological replicates. The method is based on the negative binomial distribution, with variance and mean linked by local regression. The count data were normalized on the total number of counts, taking the variance and the mean of the biological replicates into account, and subjected to further pair-wise differential expression transcriptome analysis. For each gene, we estimated the fold change from the first to the second condition (i.e. temperature growth conditions on the one hand and between species under particular growth conditions on the other hand). Significant differences in expression were identified by DESeq, and we obtained an adjusted P -value (false discovery rate <0.05; Benjamini–Hochberg adjustment) for the statistical significance of this change. For each gene, ‘normalized count data’ characterized a species for each thermal regime and have been used in some figures or for the calculation of the PI. Data are available on the Coffee Genome Hub ( http://coffee-genome.org/ ), an integrative genome information system ( Dereeper et al. 2015 ).

Statistical analyses on biomass and metabolites

The biomass and the different metabolites studied (Chl, starch, sugars, lipids and hormones) were analyzed by analyses of variance (ANOVAs) followed by the Tukey HSD (Honest Significant Difference) test at a threshold of 0.05. The correlation between metabolites and expression gene was calculated using the Person’s correlation. Correlation with a probability of P > 0.99 was considered.

Plasticity index

In order to compare the transcriptome plasticity of each species, we calculated a PI of the whole transcriptome for each TR.

$$\mathit{P.}\mathit{I}={\displaystyle \sum _{i=1}^n{\left(x_{ij}-\left(y_i\right)\right)}^2/\left({\displaystyle \sum Y_{ji}/a}\right)}$$

where x is the count data (i.e. mean of two replicates) for gene i of species j ; y the mean of counts for gene i for the three species; n = number of genes; Y = the sum of count data for species j ; a = number of species.

Gene categorization in differential expression patterns

We used a categorically partitioned analysis of the full set of genes as defined by Rapp et al. (2009) and Yoo et al. (2013) . This analysis procedure which can be used to class gene expression according to a conventional genetic model (dominance, overdominance and additivity) is usually used for a single environment. Here, we used it simultaneously for the three TRs for which we were able to compare C. arabica with its two parents (18–14, 23–19 and 28–24°C). To carry out this analysis, we first eliminated genes with low expression by fixing a limit threshold. To fix that threshold, we considered the sum of the counts (S; i.e. normalized expression levels) for the three TRs (18–14, 23–19 and 28–24°C). Only the genes for which S was ≥250 counts for at least one of the three species were kept. Indeed, beyond that threshold, we considered that we were taking a risk of detecting a signal without any biological significance. Then, we classed the genes per category for each TR. The expression level was binned into 10 patterns depending on relative expression levels between C. arabica and its parental species. The patterns were the following: ‘no change’, ‘ C. eugenioides -like dominance’ (negative or positive), ‘ C. canephora -like dominance’ (negative or positive), ‘overdominance’ (negative or positive), ‘additivity’ (negative or positive) and ‘not stable’. For a gene to be binned in the ‘no change’ category, no significant differences ( P > 0.05) had to exist between the TRs and between the species for each TR. For a gene to be classed as ‘ C. eugenioides -like dominant’ (i.e. ‘dominant C. canephora -like’, ‘overdominant’, ‘additive’ etc.), no significant difference had ever to be observed between C. eugenioides and C. arabica , but, on the other hand, there always had to be a significant difference between C. arabica and C. canephora under all three TRs and, in addition, there had to be no significant differences ( P > 0.05) between TRs. This approach covered genes for which C. arabica was in an additive, dominance or overdominance situation. The other genes were binned in the ‘not stable’ category. This represented the genes for which there could exist some significant differences between the TRs for at least one of the three species and/or the genes for which there existed significant differences between species for at least two TRs. This ‘not stable’ pattern therefore included so-called ‘plastic’ genes, while the genes of the other nine patterns were considered as ‘homeostatic’.

Funding

This work was fully supported by CIRAD.

Disclosures

The authors have no conflicts of interest to declare.

Abbreviations

Abbreviations
 
• AAO

  ascorbate oxidase

 
• ABAGE

  ABA–glucose ester conjugate

 
• APX

  ascorbate peroxidase

 
• C16:0

  palmitic acid

 
• C18:3

  linolenic acid

 
• CAB

  Chl a-b binding protein

 
• CAO

  chlorophyllide a oxygenase

 
• DPA

  dihydrophaseic acid

 
• FAD

  fatty acid desaturase

 
• FAME

  fatty acid methyl ester

 
• Fd

  ferredoxin

 
• GI

  Gigantea

 
• LHY

  Late elongated hypocotyl

 
• MAPK

  mitogen-activated protein kinase

 
• PI

  plasticity index

 
• PPO

  polyphenol oxidase

 
• PPR

  pentatricopeptide repeat

 
• PRR

  Proline-rich protein

 
• RFO

  raffinose family oligosaccharide

 
• RNAseq

  RNA sequencing

 
• ROS

  reactive oxygen species

 
• RP17

  60S ribosomal protein L7

 
• TOC1

  Timing of cab expression 1

 
• TR

  thermal regime

 
• ZT

  Zeitgeber time

 
• ZTL

  Zeitlupe

Author notes

⁷ These authors contributed equally to this work.