Temperature Impacts the Response of Coffea canephora to Decreasing Soil Water Availability

Climate change is expected to result in more frequent periods of both low rainfall and above normal temperatures for many coffee growing regions. To understand how coffee reacts to such change, we studied the physiological and gene expression responses of the clonal variety C. canephora FRT07 exposed to water deficits under two different temperature regimes. Variations in the time-dependent impact of water deficits on leaf stomatal conductance and carbon assimilation were significantly different under the 27 °C and 27 °C/42 °C conditions examined. The physiological responses 24 h after re-watering were also different for both conditions. Expression analysis of genes known to respond to water deficits indicated that drought-related signaling occurred at both temperatures. Deeper insights into the response of coffee to water deficits was obtained by RNASeq based whole transcriptome profiling of leaves from early, late, and recovery stages of the 27 °C experiment. This yielded expression data for 13,642 genes and related differential expression analysis uncovered 362 and 474 genes with increased and decreased expression, respectively, under mild water deficits, and 1627 genes and 2197 genes, respectively, under more severe water deficits. The data presented, from a single clonal coffee variety, serves as an important reference point for future comparative physiological/transcriptomic studies with clonal coffee varieties with different sensitivities to water deficits and high temperatures. Such comparative analyses will help predict how different coffee varieties respond to changing climatic conditions, and may facilitate the identification of alleles associated with high and low tolerance to water deficits, enabling faster breeding of more climate-smart coffee trees.


Introduction
Coffee is a globally important agricultural commodity that is grown mainly in developing tropical countries. World coffee production for 2017/18 is predicted to be over 9 million metric tons (USDA 2017), and much of this is produced by small holders. The revenue associated with growing green coffee supports the livelihoods of tens of millions in the tropics (Bigger and Hillocks 2007;Gennari 2015), and there is rising concern that these coffee farmers, and the economically important industry they underpin, will become increasingly impacted by climate change. Probably the most immediate consequences of climate change on coffee production will be from erratic precipitation fluctuations that are predicted to be more frequent in many coffee areas in the future, including fluctuations leading to increasing periods of prolonged water deficits (Moat et al. 2017;Grosjean et al. 2016). It is also possible that episodes of low rainfall in the coming years will also be accompanied by periods of above normal temperatures . In addition to modelling studies predicting future climate change impacts on several major coffee regions, analysis of existing weather data indicate that abnormally dry and hot periods are already occurring in some areas (Cardenas 2014;Craparo et al. 2015;Moat et al. 2017).
Given the growing recognition that climate change could significantly impact future coffee production, an increased effort has been directed towards an understanding of how, and to what extent, Arabica (Coffea arabica L.) and Robusta (Coffea canephora Pierre ex Froehner) will react to the climatic shifts that are predicted for different coffee growing areas. Much of the recent research has focused on studies concerning the response of coffee to environmental factors, including water deficits, elevated atmospheric CO 2 levels, and temperature. A smaller number of studies have also investigated indirect climate related effects on coffee, for example, how increased temperatures could lead to higher numbers of predators, like the coffee berry borer (Jaramillo et al. 2009), or the white stem borer (Kutywayo et al. 2013). Although it has been suggested that altered climatic conditions can influence the severity of coffee rust (Avelino et al. 2015), more recent work suggests the high level of coffee rust damage that occurred in Columbia between 2008 and 2011 may not be strongly linked to climate during this period (Bebber et al. 2016).
Of all the future climate-related constraints that are expected to impact coffee, the increasing severity and frequency of water deficits are predicted to have the biggest influence on production (DaMatta and Ramalho 2006; Moat et al. 2017). Consequently, many climate related experimental studies on coffee plants have focused on the effects of reduced water availability on coffee physiology, gene expression, and biochemical composition. For example, Marraccini et al. (2012) and Vieira et al. (2013) compared variations in the physiological and gene expression responses of drought tolerant and drought sensitive Robusta varieties, and suggested that the drought tolerant clones had subtle differences in ABA signaling and a better induction of the coffee antioxidant and osmoprotection systems. In another examination of the physiological effects induced by drought, (Menezes-Silva et al. 2017) explored the effect of repeated rounds of water deficit and recovery on one Robusta variety. The physiological data obtained indicated that plants exposed to multiple water deficit events had higher photosynthetic rates then control plants. They also showed that the increased photosynthetic activity was accompanied by an increase in the activity of enzymes associated with carbon and antioxidant metabolism Another form of water deficit, due to a change in atmospheric relative humidity, was recently studied by Thioune et al. (2017).
Although rising atmospheric carbon dioxide levels is a key driver of climate change, this may also increase crop plant productivity in some situations (Deryng et al. 2016). Recent work has found that higher atmospheric carbon dioxide could increase photosynthesis for coffee under normal water and nutrient availability conditions, thus suggesting future increases in global CO 2 levels may not dramatically alter coffee production Verhage et al. 2017). Other research on coffee explored the effects of both elevated CO 2 and above normal temperatures. One group of experiments indicated that higher CO 2 reduced many of the negative effects of higher temperature on photosynthesis . Related work following metabolite and enzymatic activity changes in coffee leaves at different temperatures found that adequately watered plants subjected to increased temperatures and high CO 2 levels exhibited several metabolic changes favoring acclimation to higher temperatures, although this varied between the different coffee varieties studied .
Water deficits are almost inextricably accompanied by elevated temperatures and it is well documented that both separately and in combination they can negatively impact coffee production (DaMatta and Ramalho 2006) and influence coffee quality traits (Bertrand et al. 2012). However, there are almost no studies that examine the detailed physiological and gene expression changes associated with both water deficits and higher temperatures. The work reported in this study examined the effect of progressive water deficits and recovery at 27°C and at 27/42°C on the leaf physiology of a clonal C. canephora variety, as well as providing a detailed examination of the leaf transcriptomic changes associated with water deficits in the 27°C experiment. This first RNA-Seq-based transcriptomic analysis of progressive water deficits and recovery for a clonal variety of C.canephora uncovered an extensive new set of coffee genes that exhibit differential expression in response to water availability. New research exploring these water deficit response genes, as well as related alleles of other coffee varieties, should contribute to an improved understanding of the response of coffee under natural conditions. Furthermore, as the experiments also report on recovery after exposure to water deficits, the information should also help shape ideas on the resilience of coffee after prior exposure to water deficits.

Results
In order to understand the effects of a water deficit on Robusta at both normal and elevated temperatures, we set up experiments to examine the physiological and gene expression changes in the leaves of Coffea canephora variety FRT 07 plants subjected to increased water deficits under two different temperature regimes (27°C and 27°/42°C). In the first experiment, FRT 07 plants growing at 27°C were not watered for 9 days and the response of these plants were compared to those of identical wellwatered control plants.. In the second experiment, one set of FRT 07 plants were left without watering for 9 days, but were subjected to a daily temperature regime of 27°C for 15 h and 42°C for 9 h. The corresponding controls in this second experiment were exposed to the same temperature regime but were well watered (heat treatment, no water deficit). In order to study the early part of the recovery process, plants subjected to water deficits were rewatered and the associated responses were examined.
Leaf Relative Water Content of FRT 07 plants at 27°C Versus 27°/42°C and ± Water Deficit Visual observation of the plants deprived of water at 27°C showed that the leaves of these plants looked normal at day 3 but began to show the first signs of drooping at day 7, with significant drooping seen at day 9. The leaves of plants deprived of water under the 27°C/42°C temperature regime looked normal up to day 4, with some drooping at day 6, followed by severe drooping, and some evidence of leaf desiccation, at day 9. Figure 1 shows that strong, statistically significant reductions in RWC were found for both sets of plants by day 9. Plants subjected to water deficit at the higher temperature had the greater loss of leaf water at day 9 (27°/ 42°C RWC 43.1% versus 27°C RWC 56.3%). When the water-deprived plants were re-watered, the leaf drooping associated with both the 27°C and 27°/42°C treatments was reversed completely after 4 h, and dramatic recoveries in RWC were also observed ( Fig. 1). Measurements of Fv/fm indicated that there was little change in photosynthetic light use efficiency for the plants subjected to water deficits at either 27°C or 27°C/42°C ( Supplementary Fig. 2).

Stomatal Conductance and Assimilation of FRT 07 plants with and without a Water Deficit at 27°C and 27°/42°C
As the water deficit increased with the 27°C temperature regime, g s fell between days 3 and 7, and then remained very low until day 9 ( Fig. 2A). Only a very small level of recovery in g s was observed for the plants subjected to a water deficit and grown at 27°C, 24 h after re-watering. The g s values for the well-watered control plants grown at 27°C remained similar over the 10 days of the experiment. The response of g s to water deficits at 27°/42°C was quite different (Fig. 2B). In this case, g s decreased slightly after 2 days of withholding water for both sets of plants (high temperature, with and without watering), then increased slightly at the 4 day time point, although leaf to leaf variations were very high. At day 6, g s values for both sets of plants were clearly less than those seen at days 2 and 4. Interestingly, at day 9, only the plants subjected to water deficits showed a further drop in g s , attaining a lower g s than seen in the 27°C treatment with a water deficit. It was noted that the standard deviations for both control and treated plants at 27°/42°C were generally higher than those for the plants at 27°C, suggesting more heterogeneity in the g s response at the higher temperatures. In contrast to the plants exposed to water deficits at 27°C, recovery after re-watering for the plants at 27°/42°C, with or without a water deficit, was associated with a dramatic rise in g s .
No clear effect of increasing water deficits on CO 2 assimilation (A) at 27°C was found until days 7-9, when A decreased markedly (Fig. 2C). For the control plants at 27°C the A values remained relatively stable. After re-watering, A increased significantly for the 27°C plants exposed to water deficits, although the values after 24 h of recovery were only approximately half of those seen at the start of the experiment. The response of A to water deficits for the 27°/42°C temperature regime was quite different to that seen at 27°C. Exposure to the higher temperature regime and a water deficit resulted in a dramatic reduction in A after only 2 days, and this remained low, but relatively stable, for the remainder of the period of water deficits. Interestingly the 27°/42°C control plants also showed a gradual fall in A up to day 4, after which this trend was reversed up to day 9. After re-watering and recovery for 24 h at 27°C, A increased dramatically for these 27°/42°C plants subjected to a water deficit. A small increase in A was also observed for the related control plants (high temperature treatment) after the 24-h recovery at 27°C.

Comparison of Leaf Gene Expression Changes for Plants Subjected to a Water Deficit at 27°C and 27°/42°C
The expression data obtained for previously characterized genes that are inducible by water deficits, CcDH1a and CcDREB1a, and the heat responsive gene CcHSP90-7 are shown in Fig. 3. As found previously (Simkin et al. 2008), the dehydrin gene CcDH1a transcript levels rose markedly under increasing water deficits for both temperature regimes, then the levels of CcDH1a fell rapidly and dramatically 4 and 24 h after re-watering. For the high temperature control plants, only a very weak induction of CcDH1a was detected, with the highest level of transcripts seen at day 6. CcDREB1a encodes a coffee drought inducible Dehydration-Responsive-Element-Binding transcription factor of the APETALA2 (AP2) gene family (Alves et al. 2017). The data in Fig. 3 indicates that DREB1a only showed detectable transcript increases at day 9 for the 27°C and 27°/42°C water deficit treatments, with possibly a smaller increase at day 9 for the 27°/42°C water deficit treatment. No clear increase in CcDREB1a transcripts were detected during the 27°/42°C high temperature treatment only, although some unexpected variations in CcDREB1a levels were seen at the start (T = 0) and to a lesser extent, at the end of this experiment. Examination of the transcript levels for the CcHSP 90-7 heat shock gene shows that Fig. 2 Stomatal conductance and carbon dioxide assimilation for Coffea canephora FRT07 plants subjected to a water deficit and a water deficit under a heat treatment. Squares and triangles represent control samples and water deficit samples at each temperature respectively. Panels A and B -stomatal conductance, g s , mol.m -2 .s -1 ; Panels C and D -Assimilation A, μmol CO 2 .m -2 .s -1 . Statistical analysis was performed as in the Materials and Methods. Significance levels are indicated as * P < 0.1, ** P < 0.05, ***P < 0.01 Fig. 1 Leaf Relative Water content measurements for Coffea canephora FRT07 plants subjected to a water deficit and a water deficit under a heat treatment. Two water stress experiments were carried out at 27°C and at 27°/42°C as described in the Materials and Methods. At each time point, single leaves were harvested from the center of each plant. The data presented here for both stressed and control plants represent measurements from three biological replicates ±SD. Squares and triangles represent control samples and water deficit samples at each temperature respectively. Statistical analysis was performed as in the Materials and Methods. Significance levels are indicated ** P < 0.05, ***P < 0.01 there was no significant accumulation of these transcripts during the water deficit treatment at 27°C. However, CcHSP90-7 transcript levels did increase in both the controls and water deficit treatment at 27°/42°C. In both cases, after the 24 h recovery period ended, the level of CcHASP90-7 transcripts returned to the levels seen at T = 0.

Transcriptomic Analysis of Leaves during Increasing Water Deficit and Recovery at 27°C
To identify new coffee genes exhibiting altered transcript levels under mild and severe water deficits, and to examine the response of these genes after re-watering, we employed a total RNA-Seq approach using leaf RNA associated with different stages of the 27°C water deficit experiment. The RNA-Seq libraries generated, which were made using ribosomal RNA-depleted total RNA samples, each yielded over 70 million reads (Table 1). After quality checks, the read data from the 4 time-points of this experiment were mapped directly onto the annotated Robusta genome (Denoeud et al. 2014), resulting in relative expression values (CPM -counts per million) for 13,642 genes.

Identification of Genes with Strong Differential Expression under Mild and Strong Water Deficits at 27°C
Differential gene expression analysis was carried out to identify the FRT 07 genes with significant transcript level changes as the water deficit increased, that is, transcript levels at T = 0 versus those found at T = 1 (3 days, mild water deficit) and T = 3 (9 days -stronger water deficit). The Venn diagram presented in Fig. 4 shows that transcript levels increased or decreased more than 2-fold (log 2 1; P < 0.05) for 362 genes and 447 genes, respectively, in the leaves under the mild water deficits (T = 1). As expected, the stronger water deficit at T = 3 caused more variations in gene expression, with transcript levels increasing for 1627 genes and decreasing for 2197 genes, respectively. Fig. 4 also shows that 123 genes were uniquely up regulated under mild water deficits, while 892 genes were uniquely up regulated as a result of the stronger water deficit, with 191 genes up regulated under both conditions. Similarly, 135 genes were uniquely down regulated at T = 1, while 1423 genes were uniquely down regulated at T = 3. There were 283 genes down regulated under both levels of water deficit. Quantitative expression analysis of genes CcDH1a, CcDREB1a, and CcHSP90-7 in leaves during the different water deficit and/or heat treatments. Quantitative RT-PCR gene expression profiles were determined for RNA extracted from one-half of the leaves used for the RWC analysis in Fig. 1 as described in the Materials and Methods. Transcript levels (n = 3 ± SD) were determined relative to the expression of the gene ribosomal protein like 39 (rpl-39) in the same sample. The plants subjected to water stress were re-watered after the 9 day sample was taken. Statistical significance was determined for each time sample versus the T = 0 sample as described in the Material and Methods. * P < 0.1, ** P < 0.05, ***P < 0.01 The RNA-Seq data not only allowed us to identify genes whose expression changes in leaves during a water deficit at 27°C, it also enabled us to follow the changes in transcript levels for each of these genes across two levels of water availability, as well as after re-watering (T = 5, 24 h after rewatering). To examine these changes in more detail, we examined some of the genes that exhibited the biggest UP changes in transcript levels at T3. The genes chosen are listed in Supplementary Table 2, and include several transcription factors (TFs) such as NAC (2), MYB39 related (2), ethylene responsive transcription factor (1), heat stress related TFs (2), WRKY (8), and an NFYA type transcription factor, which exhibited a quite large transcript increase. We also found three DREB related genes with increased transcripts at T3. The sequence of one of these T = 3 induced genes (DREB-1D (Alves et al. 2017)) is nearly identical to CcDREB1a that we have shown earlier to be induced at T = 3 by qRT-PCR ( Fig.  3). We also found several "transport protein" associated genes with increased transcript levels at T = 3, for example one Na+/ Ca2+ exchanger, a boron transporter, an aluminum activated malate transporter 2, a monosaccharide-sensing protein 2, and a sugar transport protein. A number of hydrolase annotated genes similarly showed increased transcript levels at T = 3, such as two endo-beta-glucanases with relatively large transcript increases, as well two xyloglucan endotransglucosylase/ hydrolases, and a neutral invertase.
Other genes with notable increases in transcript levels at T = 3 included an asparagine synthetase [glutamine-hydrolyzing], and several ethylene synthesis related genes (two 1aminocyclopropane-1-carboxylate synthases, and several 1aminocyclopropane-1-carboxylate oxidases). We also found significant transcript increases for defense related genes like momilactone A synthase, a gene encoding a phytosulfokine precursor (a possible defense related signaling molecule), a gene for an aldehyde dehydrogenase, a gene associated with late blight resistance, and a gene encoding a member of the NDR1/HIN1-like gene family. A few genes associated with phenolic metabolism also exhibited significant transcript level increases at T = 3, such as two members of the 4-coumarate-CoA ligase gene family, two genes for bifunctional 3dehydroquinate dehydratase/shikimate dehydrogenase, which reversibly converts quinate to shikimate, and two anthocyanin 3'-O-beta-glucosyltransferases. Interestingly, several ABA associated genes showed increased transcript levels at T = 3, such as those that encode an ABA degrading enzyme (abscisic acid 8′-hydroxylase 4), an MLO-like protein 6, an NDR1/ HIN1-like 2 protein, and two genes encoding "aspartic proteinase in guard cell-1" protein.
Examples of genes showing strong reductions in transcript levels at T = 3 (Supplementary Information Table 3) include multiple MYB and bHLH related transcription factor genes. A decrease in transcript levels were also observed for metabolite transport genes, such as several encoding ABC transporter proteins, three members of the cationic amino acid transporter family, a metalnicotianmine transporter gene, a major facilitator superfamily protein (putative auxin transporter), and several encoding MATE efflux family proteins. A number of hydrolase related genes also had strong reductions in transcript levels at T = 3, such as two genes encoding polygalacturonase-1 non catalytic subunit beta (one with very high expression), two pectate lyases, and three chitinase related genes. A few genes encoding enzymes potentially involved in the metabolism of phenolic compounds also had reduced transcript levels at T = 3, again suggesting a reorganization of phenolic metabolism during episodic water deficits. For example, transcript reductions were seen for three HXXD/BAHD type acyltransferases, three bifunctional dihydroflavonol 4 reductases, and two genes annotated as caffeic acid 3-O-methyltransferases. Other genes exhibiting striking reductions in transcript levels at T = 3 include a group linked to chlorophyll synthesis, including several highly expressed genes encoding chlorophyll a-b binding proteins and two magnesiumchelatase subunit H genes. Transcript reductions were also seen for multiple genes containing BURP domains, including a number of dehydration responsive protein genes (RD22). Marked transcript reductions at T3 were also found for genes annotated as "auxin related" and "aspartic endoproteinase in guard cells". Finally, significant expression decreases at T3 were observed for two ABA receptors (PYL6 and PYL1), three aquaporins, four HSP genes, a carbonic anhydrase, and a putative carotenoid cleavage dioxygenase (both the latter two genes were highly expressed in the leaves of well-watered plants. Figure 4 indicates that 388 genes and 290 genes showed unique increases or decreases of transcript levels >2 fold (log 2 1; P < 0.05) after 1 day of recovery from water deficits (T = 5). Of the genes showing transcript increases uniquely at this stage, only 23 of these showed transcript increases >8 fold (log 2 3) in comparison to T = 0 (Supplementary Information Table 4). Five of these genes were unknown genes (including three of the top four highly expressed genes of this group). This group included genes encoding two transcription factors (a NAC and an ethylene responsive factor), two LRR receptor like serine/threonine protein kinases, and four genes potentially involved in cell wall modifications (two xyloglucan endotransglucosylase/hydrolases, a pectin-methylesterase/ PME inhibitor, and a cationic peroxidase). Other genes with a uniquely high expression at T5 were two cytochrome p450 genes, a phenylalanine ammonia-lyase (PAL) gene, a salicylate O-methyltransferase, a BON association protein, and a 1deoxy-D-xylulose-5-phosphate synthase potentially involved in isoprenoid synthesis. Of the genes showing transcript decreases at T5, only three exhibited >8 fold lower expression in comparison to T = 0 ( Supplementary Information Table 4), and two of these genes were unknown genes, and the third encoded histone H2B.2. Only nine additional genes showed >5 fold lower expression at T5. Five of these genes encode uncharacterized proteins, the others are a monothiol glutaredoxin, a BTB/POZ/TAZ domain encoding protein, an IQ domain containing protein, and a Tesmin/TOS1-like CXC domain containing protein.

Gene Onthology Enrichment Analysis of 27°C RNA-Seq Data for Mild and Strong Water Deficits and for 24 Hours after Rewatering
Heatmaps for the UP and DOWN differentially expressed gene sets were generated from a "Biological Processes" Gene Onthology (GO) based enrichment analysis. Fig. 5A shows some of the processes "UP only at T=1" include two processes each for ABA metabolism, alcohol metabolism, inositol metabolism, NAD metabolism, and organic hydroxy compound metabolism. Other uniquely "UP only at T= 1" genes were for polysaccharide catabolism, proline catabolism, and for reactive nitrogen species metabolism, as well as for responses to oxidative stress and sulfate reduction. Processes that are "UP only at T=3" include ceramide metabolism, clatherin assembly, dephosphorylation, and a glutamine family amino acid related process. Other "UP only at T = 3" genes were for glycoprotein and glycosphingolipid metabolism, nucleobase containing metabolism, organic cyclic biosynthesis, two for phospohtidylserine metabolism, three for response processes (to oxygen containing compounds, to water, to water deprivation), and two for signaling related processes. Finally, "UP only at T=5" processes include two each for genes encoding aromatic amino acid metabolism, isoprenoid metabolism, L-phenylalanine metabolism, multiorganism/multicellular related, nicotianamine, and reproductive processes. Other "UP only at T5" processes found were for erythrose 4 phosphate/phosphoenolpyruvate, exocytosis, external encapsulating structure organization, organic acid catabolism phosphorylation, protein phosphorylation, TCA biosynthesis and urea related metabolism.
Some interesting "DOWN only at T=1" processes ( Fig.  5B) include two for multi-organism/multicellular organisms, nicotianamine metabolism, transport (organic acids/ organic anions), and reproduction. Other "DOWN only at T = 1" included asparagine biosynthesis, cell recognition, drug membrane transport, glutamate metabolism, isoprenoid biosynthesis, lipid metabolism, oxidation-reduction process, response to oxidative stress, transition metal ion transport, TCA biosynthesis, and vitamin biosynthesis. Notable "DOWN only at T=3" processes were two for the mitotic cell cycle, protein polymerization/refolding, as well as multiple cellular metabolism related processes and glycogen related processes. The final group of "DOWN only at T = 3" processes were related to amino sugar biosynthesis, biological adhesion, cell division, cell septum assembly, cellular component biogenesis, hormone metabolism, regulation of biological quality, and regulation of hormone levels. "DOWN only at T=5" (recovery only) processes included two each for chlorophyll biosynthesis/metabolism and monocarboxylic acid biosynthesis/metabolism, and one process each for base excision repair, cellular polysaccharide biosynthesis, cellular response to stimulus, ion transport, and RNA processing.

Discussion
A series of climate modelling studies have predicted that climate change will increasingly threaten coffee production world-wide due to more frequent occurrences of above average temperatures and longer periods with little or no rainfall (Läderach et al. 2017;Bunn et al. 2015;Ovalle-Rivera et al. 2015). Whilst there are several studies on the effects of water deficits and high temperatures on coffee there are no studies, to date, that have focused on the combined effects of both water deficits and higher temperatures, although these are climate related constraints that often occur together. We addressed this knowledge gap by examining the physiological and gene expression responses of a clonally propagated Robusta variety to two different temperature regimes (27°C versus 27°C/42°C), with or without a water deficit. Examination of the physiological responses for the 27°C treatment with a water deficit showed a clear decrease in RWC by day 9, relative to the well-watered control plants; followed by a return to near control RWC levels 4 h after rewatering. The water deficit at 27°C was also associated with a dramatic reduction in stomatal conductance (g s ) by day 7 compared to the controls, and this value remained low to day 9. Unexpectedly, little recovery in g s was observed for the 27°C water deficit treated plants 24 h after re-watering, indicating some residual impact on stomatal movements and their sensitivity to increased water availability. In grapevine, a similar failure of stomata to recover rapidly after a prolonged water deficit has been attributed to foliar accumulation of ABA, and it was suggested that the slow recovery in g s may favor the repair of xylem embolism damage (Tombesi et al. 2015). For plants subjected to a water deficit at 27°C A remained similar to the watered controls up to day 7, after which it declined. Upon re-watering, A recovered partially, but this occurred in the absence of a significant recovery in Fig. 5 GO enrichment heatmap of UP and DOWN regulated genes at different stages of the 27°C water deficit treatment. Differentially expressed "UP" and Down genes were subject to enrichment for GO biological processes. Enriched processes with an FDR <0.05 were included in the heatmap. Panel A -"UP" Regulated genes, Panel B "Down" Regulated genes. Expression changes are relative to T = 0 and samples presented correspond to T = 1 (3 days), T = 3 (9 days), and T = 5 (24 h after water added -recovery) g s . This suggests that the recovery of A after a water deficit may be due, in part, to diffusional modifications in CO 2 transport that are independent of g s (Urban et al. 2017).
During the 27°/42°C high temperature only experiment, there was no drop in leaf RWC, but a significant drop in RWC was seen during the combined high temperature and water deficit experiment between days 6 and 9. Interestingly, g s for both the high temperature alone and the high temperature with water deficits treatments exhibited a similar downward trend for days 2-6, after which time g s stabilized for the plants at high temperature, but fell further for the plants subjected to the combined high temperature and water deficit treatment (day 9). Both sets of plants at 27°C/42°C showed a fast return to the initial g s values after the 24 h recovery period at 27°C, contrasting with the poor recovery of g s seen for the plants subjected to a water deficit at 27°C, suggesting temperaturerelated modifications in ABA production or re-mobilization (Tombesi et al. 2015). Alternatively, the higher temperature may have altered ABA and/or other signaling pathways in the leaf in ways that facilitate a rapid recovery in g s after the temperature was reduced and the water deficit removed.
Values for A of plants subjected to high temperatures alone showed a clear downward trend at days 2-4, but then increased from days 6-9, indicating an acclimatization had occurred without any alteration in g s . A further rise in A occurred during the 24 h recovery period at 27°C, and this increase in A was accompanied by a rise in g s . Plants subjected to the 27°/ 42°C heat treatment and a water deficit exhibited a larger decrease in A at day 2 than seen for the 27°/42°C controls, and A remained low up to day 9. After re-watering, the A values for these plants increased dramatically as did g s .
Overall, it appears that from day 4, the photosynthetic metabolism of the plants exposed to the higher temperature appeared to show some adaptation in A, but this capability was lost when combined with a concomitant water deficit. Whilst an ability to acclimate to different environmental constraints is a common feature of a number of plant species (Lamke and Baurle 2017), such adaptions may be limited under field conditions with the addition of further environmental constraints. The rapid recovery of A after re-watering for both sets of 27°C/42°C plants indicates that many of the effects of high temperature and water deficit are reversible, although the recovery may depend on the duration and frequency of the exposure. Higher assimilation levels in plants subjected to repeated cycles of varying water availability (Menezes-Silva et al. 2017) suggest that coffee may even exhibit improvements in photosynthesis after recurrent water deficits.
It has been shown that tomato varieties with differences in heat tolerance exhibit small, but specific physiological response variations to heat, drought and both stresses together (Zhou et al. 2017). Therefore, we decided to examine whether a second coffee variety (FRT 06) might show different timedependent trends in g s and/or A in response to falling soil water availability at different temperatures than FRT 07 using the methodology presented here. This preliminary experiment indicated that FRT 06 plants exhibited mostly similar timedependent response trends for RWC, g s , and A versus FRT 07 with and without water (at 27°C and 27°/42°C for 9 days; E.H. Thioune, unpublished results). However, we did detect two response differences. First, the A values for the wellwatered FRT 06 plants at 27°/42°C fell similarly to the FRT 06 plants at 27°/42°C subjected to a soil water deficit. Second, the well-watered FRT 06 plants at 27°/42°C appeared to lack the surprisingly strong adaptation of A seen for the well-watered FRT 07 plants after day 4 under the 27°/42°C (heat only; Fig. 2D). Such observations suggest future work comparing time-dependent RWC, g s , and A responses of clonal Robusta varieties with different water deficit and/or heat responses may enable the identification of coffees with contrasting physiological responses to these stresses.
To obtain transcriptomic data for the water deficit experiment at 27°C, we carried out a deep RNA-Seq experiment for leaf samples taken at three of the five time sampling periods (T = 1, T = 3, T = 5). Analysis of the RNA-Seq data uncovered a large number of differentially expressed genes (DEGs versus T = 0) whose transcript variations can be associated with physiological changes in response to increasing water deficits. As there are too many genes to discuss in detail here, we have focused on some of the most interesting DEGs found, many of which showed >3 fold changes (>3 FC; >log 2 3) UP or DOWN at day 9 when the plants were subjected to a strong water deficit (Supplementary Tables 2 and 3). The first notable set of genes found were members of the DREB transcription factor family (Lata and Prasad 2011). Two of these had significant UP expression at day 9 (Cc8_g13960, FC = 7.85 and Cc2_g03430, FC = 5.51),with DNA sequence analysis showing the latter gene is the well characterized drought stress regulated gene CcDREB-1D (Alves et al. 2017). Recent detailed studies of the CcDREB-1D gene in different varieties found variety-specific promotor variations linked to phenotypic differences in drought tolerance , and further work found that haplotype specific CcDREB-1D promoters can influence the type and level of a reporter gene's response to water and temperature-related stresses (Alves et al. 2017;de Aquino et al. 2018). Although both DREBs noted here had very low transcript levels without a water deficit, and high levels under a strong water deficit (day 9), they appeared to react differently after re-watering, with the transcript levels falling faster for Cc8_g13960 versus Cc2_g03430 (DREB-1D), suggesting these two genes may have some differences in function. A third Up regulated DREB exhibited a relatively high expression at T = 0 (no water deficit), but had only a slight increase in transcripts at day 9 (Cc08_g09520, FC = 1.58) and a return to near T = 0 levels after re-watering. In contrast, the fourth differentially expressed DREB had a relatively significant DOWN regulation at T = 3 (Cc06_g10260, FC = −3.32), followed by a restoration to T = 0 transcript levels after re-watering.
Another interesting group of DEGs were ABA/drought signaling genes like Cc02_g18810 (FC = 2.95) whose transcript levels rose at day 9, then fell after rewatering to levels below those seen at T = 0. Cc02_g18810 is annotated as an ABA degrading cytochrome p450 type enzyme called abscisic acid 8′-hydroxylase involved in ABA catabolism and the regulation of ABA levels both during water deficits and during recovery (Liu et al. 2014). The ABA related gene Cc02_g11510, which encodes an MLO-like protein, also showed a large increase in transcript levels at day 9 (FC = 5.61), then dropped to T = 0 levels after re-watering. Lim and Lee (2014) have described an MLO-like gene in pepper that is induced by both abiotic and biotic stress, and indicated that this gene may act as a negative regulator of ABA signaling that results in a decrease in water loss from leaves under water deficits (Lim and Lee 2014). Transcriptional alterations were also seen for several putatively ABA linked NDR1/ HIN1-like (NHL) genes. One, the intron-less gene Cc06_g06450 showed significant UP expression at day 9 (FC = 4.76), followed by a drop in expression after recovery. An Arabidopsis NDR1 gene induced by ABA has been reported to be involved in the osmotic stress response (Bao et al. 2016). Intriguingly, another intron-less NDR1/HIN1 like gene (Cc06_g06460) is located on chromosome 6 right next to Cc06_g06450, but this second gene, with clearly a different sequence, showed no significant transcription differences during the experiment. Five differentially regulated ABA related genes called "aspartic proteinase in guard cell-1" were also found. Two other interesting ABA signaling related DEGs are two carotenoid dioxygenases (Cc05_g10210, FC = 1.66; Cc08_g05610, FC = −3.36). The gene Cc05_g10210 encodes CcNCED3, a gene involved in ABA synthesis (Kalladan et al. 2019) and previously shown to have increased transcript levels in coffee leaves under different water deficits (Simkin et al. 2008;Thioune et al. 2017). Transcript levels for carotenoid dioxygenase Cc08_g05610, were relatively high initially, dropped significantly as the water deficit increased, then exhibited only a partial recovery after re-watering. A third carotenoid cleavage gene DEG (Cc02_g28080) was found, but its transcript levels were not affected by water deficits. Further interesting ABA related DEGs included a PYL/ABA receptor gene (Cc02_g05990) that was significantly DOWN regulated (FC = −4.12), and returned to values above those found originally after recovery. Four other PYL receptor family genes were also noted, but only one had any major alteration in transcript levels at the higher water deficits (Cc08_g02750, FC = -1.76; Cc08_g15960, FC = 0.05; Cc02_g39180, FC = 0.59; Cc02_g01800, FC = -0.08).
An annexin with a high T = 0 expression showed a strong transcript level rise at the higher water deficit (Cc10_g07300, FC = 2.5). Annexins have been linked to drought and oxidative responses, potentially via ABA signaling (Ijaz et al. 2017) and this coffee gene was most homologous to Annexin 1 of Arabidopsis, which is a calcium transporter gene known to be induced by water deficits and high temperatures (Huh et al. 2010;Liao et al. 2017). The two other annexin DEGs (Cc07_g15280 and Cc11_g05800) exhibited only slight transcript increases. Gene onthology (GO) enrichment analysis confirmed that genes associated with the Biological Process "response to oxidative stress" were UP during low and high water stress, while others in this process were also DOWN under the high water deficit. At the level of individual genes, several potentially oxidation stress related peroxidases exhibited significant transcript level increases (>FC = 2 to FC = 6) under the high water deficit. Also, two superoxidase dismutase annotated genes showed transcript decreases of > FC = − 2 at the higher water deficit condition tested. These observations further indicate significant transcriptional changes occurred for genes involved in addressing alterations in the levels of reactive oxygen species associated with water deficits (Qi et al. 2018). Another consequence of falling soil water availability is the development of embolism in xylem vessels (Klein et al. 2018). As embolism may have begun to occur in the leaf xylem as the water deficit progressed, we looked for significant expression changes of genes associated with embolism, such as sucrose transporters and aquaporins (Secchi and Zwieniecki 2011). Interestingly, as the water deficit increased, a very strong increase in expression was observed for one sugar transporter (Cc02_g32480, FC 5.02), and significantly decreased expression was seen for one PIP aquaponin gene (Cc02_g16640, FC -3.77) and two TIP aquaporins (Cc08_g09400, FC -4.081 and Cc01_g12260, FC -1.45). Future experiments, such as looking at the localized expression of these genes in the xylem vessels of coffee leaves and stems during water deficits should be informative.
The phytohormone ethylene is involved in many plant signaling pathways, including those related to abiotic and biotic stresses (Dubois et al. 2018). We found a large UP transcript i o n a l c h a n g e f o r t w o g e n e s a n n o t a t e d a s 1aminocyclopropane-1-carboxylate synthases (ACS), which encode the enzyme making the precursor of ethylene 1aminocyclopropane-1-carboxylic acid (ACC) and conjugates of ACC (Van de Poel and Van Der Straeten 2014). Transcript levels for these ACS genes (Cc07_g07490, FC = 10.84 and Cc08_g15220, FC = 4.42) were nearly undetectable initially but rose as the water deficit increased and then fell back after rewatering, although the transcript levels of both genes remained well above the initial levels. Several DEGs found were annotated as ethylene synthesis enzyme 1aminocyclopropane-1-carboxylate oxidase (ACO), and some showed relatively significant UP expression at day 9 (Cc11_g07020, FC = 2.58; Cc05_g02920, FC = 2.21; Cc05_g02900, FC = 1.85; Cc05_g02890, FC = 1.56). We note that the transcript levels of Cc11_g07020 returned to near original levels after re-watering, but the levels of the other three, which are clustered together on chromosome 5, actually increased further after re-watering. Another ACO gene showed relatively high expression prior to the water deficit, but declined as the water deficit increased (Cc09_g01110, FC = −2.24 day 9). We also found other ethylene signaling genes with significantly altered expression at day 9 (T = 3), including a predicted ethylene dependent gravitopismdeficient like gene (Cc02_g06740; FC = −3.13), multiple ethylene-responsive transcription factors (Cc05_g07590, FC = 2.21; Cc02_g03420, FC = 6.98; and Cc06_g12520, FC = 2.96), and two AP2 related TF genes (Cc10_g07680; FC = −3.09 and Cc08_g15690, FC = −3.13). Transcript changes for some ethylene related TFs could be due to overlaps between wound/damage-related responses and drought related responses (Heyman et al. 2018). This is supported by the observation that several defense related genes also showed transcriptional changes at the higher water deficit, including putative disease resistance proteins (Cc01_g02600, FC = −5.16 and Cc02_g26510, F = 3.95) and pathogen related proteins (Cc11_g07780, FC = 4.68 and Cc10_g14270, FC = 4.3). A search for DEGs annotated as ethylene receptors uncovered four genes, but only one (EIN4) showed significantly altered expression at day 9 (Cc03_g05070, FC = −1.37). It has been noted recently (Bakshi et al. 2018) that the Arabidopsis EIN4 is involved in "crosstalk" with ABA signaling.
Water deficits are known to induce genes associated with chlorophyll degradation and metabolite recycling (Zhu et al. 2017), and these processes eventually lead to cell death. Thus it was not surprising we found large transcript decreases for genes associated with chlorophyll reduction, such as decreased transcript levels for several highly expressed light-harvesting chlorophyll a/b-binding (LHCB) proteins like Cc04_g16410 (FC = −5.52) and a gene encoding a putative magnesium-chelatase subunit H protein (Cc06_g17100, FC = −4.15). A particularly large increase in transcript levels was observed for an asparagine synthetase gene (Cc01_g14420, FC = 6.76), suggesting that a switch to nitrogen recovery/recycling may have been triggered to recover nitrogen from the dehydrating leaves (Masclaux-Daubresse et al. 2010;Canales et al. 2012). Consistent with the idea that early senescence signaling has been activated by a strong water deficit, we also found a very large increase in transcripts for the gene NYE1 (Cc08_g13770, FC = 6.01), whose protein is directly involved in chlorophyll degradation (Ren et al. 2007;Jibran et al. 2015). We previously did not detect NYE1 induction during a "humidity shock" at 23°C when the leaves experienced significant dehydration (Thioune et al. 2017), suggesting that there may be distinct signaling differences between a water deficit associated with a reduction in atmospheric humidity and that caused by a reduction in soil water availability.
After re-watering there were unique increases or decreases in gene transcript levels associated with the 388 and 290 Robusta varieties (Fig. 4). However, relatively few of these "recovery specific" genes had transcript level increases >3 FC in comparison to the original levels. Of the 23 such genes found, three of the top four highly expressed genes were unknown genes ( Supplementary Information Table 4). Similarly, of the unique genes with lowered expression after re-watering, only three exhibited < −3 FC in expression versus the transcript levels at the start of the experiment. It will be interesting in the future to explore whether there are orthologous genes in other plants that show a similar transcription response during the recovery from water deficits.
The RNA-Seq based transcriptional data for Robusta FRT 07 identified a large number of new water deficit responserelated genes, in addition to giving more details about the responses of several previously characterized water deficit inducible genes in coffee. Furthermore, the work provides the basis for the establishment of a linkage between leaf physiological changes induced by mild and strong water deficits and the associated global changes in gene expression using a deep RNA-Seq transcriptomics approach. Future work on the newly identified water deficit response-related genes found in this work, both annotated and unannotated, should contribute to a better understanding of the impacts of water deficits on coffee, and the recovery after a water deficit. It could be argued that the rapid recovery from a water deficit may be the more important property given that some reduction in performance, at even moderate reductions in soil water, is largely unavoidable. Finally, similar studies using specific clonal varieties of Robusta with strong variations in their responses to water deficits should enable the establishment of linkages between allelic variations and transcriptomic expression for the new coffee water deficit response genes identified here and the phenotypic data, including physiological response data for these clonal varieties. Such an approach will facilitate the identification of alleles and haplotypes that either positively, or negatively, impact on the control of a plant's water balance, providing breeders of both Robusta and Arabica coffee with information to accelerate the breeding of more Climate Change resilient coffee varieties.

Plant Material
Coffea canephora Pierre ex A. Froehner (Robusta coffee) genotype FRT07 was grown from somatic embryos produced in bioreactors and the plants were grown in 0.5 L pots in a greenhouse maintained at 25-30°C with only natural daylight as described in detail previously (Thioune et al. 2017). To ensure the plants were not subjected to occasional water deficits, the plants were additionally placed inside a plastic tunnel and regularly watered (constant humidity of >80-85%). The plants used were approximately 8-16 months old with an average of at least 8-12 well developed leaf pairs.

Water Deficit Treatments
Six greenhouse grown plants were transferred to a growth chamber set up at 27°C, RH 60%, and a daily regime of 16 h light (50 μmol.m −2 .s −1 ) and 8 h dark. To acclimatize the plants to this lower humidity level, they were watered copiously under these conditions for 9-10 days prior to starting the experiments. The effects of decreasing soil water availability was followed under two different temperature conditions (each time, 3 plants subjected to water deficit and 3 control plants were continuously watered). Under condition 1, the growth chamber was maintained at 27°C and RH 60% over the course the experiment, while under condition 2, the plants were subjected to a "daytime" high temperature stress, (caused by raising the temperature to 42°C 1 h after initiation of the light phase each day), then returning the temperature to 27°C after 9 h. The RH was also maintained at 60% over the course this latter experiment. After the re-watering these high temperature stressed plants at day 9, the chamber was left at 27°C for the full 24 h recovery period. The overall experimental plan is outlined in Supplementary Fig. 1.

Physiological Measurements
At each time point, the physiological measurements were made on each plant in the incubation chamber, starting in the same order at same time each morning (approximately 1 h into the light cycle), first on the drought stressed plants, then moving on the watered control plants. The time points for the 27°C experiment were T = 0, T = 3 days, T = 7 days, T = 9 days, T = 9 days +4 h (15 min after rewatering; the stressed set only), and T = 10 days. The time points for the 27/42°C experiment were T = 0, T = 2 days, T = 4 days, T = 6 days, T = 9 days, T = 9 days +4 h (stressed plants only), and T = 10 days. The measurements were taken for one tagged leaf of each plant, which was always from the second significant pair from the apex. The fluorescence measurements were done first, using a portable chlorophyll fluorimeter (HandyPea, Hansatech, King's Lynn Norfolk, UK). The leaves were dark adapted for 10 min, using the leaf-clips supplied by the manufacturer and the F v /fm values were measured after applying a s h o r t p u l s e o f s a t u r a t i n g l i g h t ( o n e s e c o n d a t 1500 μmol.m −2 .s −1 ). This was followed by measurements of stomatal conductance (g s ) and CO 2 assimilation (A) measurements on the same leaf using a Li-COR 64000XT (Li-Cor, Lincoln, NE, USA) infra-red gas analyzer with a saturating p h o t o s y n t h e t i c a l l y a c t i v e r a d i a t i o n ( PA R ) o f 1000 μmol.m −2 .s −1 , an ambient CO 2 concentration of 400 μmol.m −2 .s −1 , a chamber relative humidity of 30-40% and a leaf temperature of 27°C. Measurements of g s and A were recorded when the readings reached stable levels.

Measurements of Relative Water Content (RWC) and Collection of Tissue for RNA Extraction
After the physiological measurements had been made for each plant, a leaf from that plant was harvested from a roughly midway position. Each leaf was divided in half along its midrib axis using a scalpel and one half was placed in a 25 ml falcon tube, which was stored on ice before determining the relative water content (RWC) as described earlier (Thioune et al. 2017). The other half was put into another falcon tube and immediately frozen in liquid nitrogen, followed by storage at -80°C for later RNA extraction.

RNA Extraction
Leaf samples were ground into a powder using liquid nitrogen and total RNA extracted using the RNeasy Plant Mini Kit (Qiagen®) following the suppliers protocol, with the addition of an on-column RNase-free Dnase I treatment (Qiagen®). Quality and quantity assessment of the total RNA was performed using the Agilent 2100 Bioanalyzer®.

Quantitative RT-PCR
The method was as described previously (Thioune et al. 2017). Briefly, cDNA was made with the Superscript III Reverse Transcriptase kit (Invitrogen, Carlsbad, California, US) as follows: 0.5 μg total RNA sample was added to 1 μl oligo dT(18) (Sigma-Aldrich, St. Louis, Missouri, US) and made up to a final volume of 12 μl with sterile RNase-free water, then incubated at 65°C for 10 min and rapidly cooled on ice. Subsequently, 4 μl of 5× first strand buffer, 1 μl of 0.1 M DTT, 1 μl of RNase Out Ribonuclease Inhibitor (Invitrogen), 1 μl of dNTP mix (10 mM each) and 1 μl SuperScript III Rnase H-Reverse transcriptase (200 U μl − 1) were added to make up a final volume of 20 μl. The reactions were incubated at 50°C for 60 min, followed by an enzyme inactivation step by heating at 70°C for 15 min, followed by storage at −20°C. The quantitative reverse transcriptase-PCR (qRT-PCR) reactions were carried out after diluting the cDNA 50 fold using TaqMan probes as described previously (Lepelley et al. 2007). Sequences of the primers and probes are given in SI Table 1. Analysis of the Ct values of the reference gene CcRPL39 for all the samples indicated that the expression of this gene was not significantly influenced either within a single treatment, or between the different treatments.

RNA-Seq Analysis
RNA sequencing was performed using the same total RNA samples as for the qRT-PCR experiments. As well as assessing RNA quality with an Agilent 2100 Bioanalyzer®, quality was further analyzed using a Fragment Analyzer (Advances Analytical) and the RNA quality was found to be homogeneous within samples. Ribosomal RNA depletion was then performed on the samples using the Ribo-Zero magnetic kit (Illumina) for plant leaves and RNA-Seq libraries were prepared from the rRNA depleted samples using the TruSeq Stranded Ribo-Zero Gold Kit (Illumina), followed by a PCR amplification step. The procedure was automated on a Sciclone NGS Workstation (Perkin Elmer). The manufacturer's protocol was followed, except for the PCR amplification step. The latter was run for 13 cycles using the KAPA HiFi HotStart ReadyMix (Kapa BioSystems). This optimal PCR cycle number has been evaluated using the Cycler Correction Factor method as previously described (Atger et al. 2015). Purified libraries were quantified with Picogreen (Life Technologies) and the size pattern was controlled using the DNA High Sensitivity Reagent kit on a LabChip GX (Perkin Elmer). Libraries were then pooled and clustered at a concentration of 8 pmol on a high output paired end sequencing flow cell (Illumina). Sequencing was performed for 2 × 125 cycles on a HiSeq 2500 (Illumina) with V4 chemistry strictly following Illumina's recommendations. Demultiplexing of the sequence data was done using Cassava (http://cancan.cshl.edu/labmembers/gordon/fastq_illumina_ filter/#contact). Ribosomal RNA was filtered using sortmerna v 2.1b (Kopylova et al. 2012). Reads were cleaned with trimmomatic v 0.36 (Bolger et al. 2014) with the following parameters: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:15 MINLEN:50. Cleaned reads were mapped to the indexed reference genome sequence of Coffea canephora (Denoeud et al. 2014) using hisat2 (Kim et al. 2015). A counts table was generated from the mapped read files using stringtie (Pertea et al. 2015) and prepDE.py script with the reference genome annotation file. Differential gene expression analysis was performed using an FDR of <0.05 with edgeR (Robinson et al. 2010) using an R script (https://github.com/bcbc-group/manuscripts/blob/master/ McCarthy_2018/Jim_edge.R). All time points were compared to T = 0. Venn diagrams were created (http:// bioinformatics.psb.ugent.be/webtools/Venn/) to illustrate overlapping sets of filtered, differentially expressed genes between T = 0 and the different time points. Gene ontology analysis was done using topGO (Alexa A, Rahnenfuhrer J (2016). topGO: Enrichment Analysis for Gene Ontology. R package version 2.32.0.) implemented through an R script (https://github.com/bcbc-group/ manuscripts/blob/master/McCarthy_2018/Jim_GO.R).

Statistical Analysis
Statistical analyses of the data was carried out using SigmaPlotR 12.0 (Systat, San Jose, CA, USA). A one-way analysis of variance (ANOVA) was applied, followed by a post hoc 'Tukey test' in order to obtain pairwise multiple comparisons among data at the same time points but at different temperatures (physiological measurements), or, between T = 0 and the other time points for each gene at each temperature (gene expression data).