Overexpression of Water-Responsive Genes Promoted by Elevated CO2 Reduces ROS and Enhances Drought Tolerance in Coffea Species

Drought is a major constraint to plant growth and productivity worldwide and will aggravate as water availability becomes scarcer. Although elevated air [CO2] might mitigate some of these effects in plants, the mechanisms underlying the involved responses are poorly understood in woody economically important crops such as Coffea. This study analyzed transcriptome changes in Coffea canephora cv. CL153 and C. arabica cv. Icatu exposed to moderate (MWD) or severe water deficits (SWD) and grown under ambient (aCO2) or elevated (eCO2) air [CO2]. We found that changes in expression levels and regulatory pathways were barely affected by MWD, while the SWD condition led to a down-regulation of most differentially expressed genes (DEGs). eCO2 attenuated the impacts of drought in the transcripts of both genotypes but mostly in Icatu, in agreement with physiological and metabolic studies. A predominance of protective and reactive oxygen species (ROS)-scavenging-related genes, directly or indirectly associated with ABA signaling pathways, was found in Coffea responses, including genes involved in water deprivation and desiccation, such as protein phosphatases in Icatu, and aspartic proteases and dehydrins in CL153, whose expression was validated by qRT-PCR. The existence of a complex post-transcriptional regulatory mechanism appears to occur in Coffea explaining some apparent discrepancies between transcriptomic, proteomic, and physiological data in these genotypes.


Introduction
Drought events have become more frequent, severe, and erratic nowadays, affecting the quality and yield of most crops [1,2]. Under the initial stages of drought, stomatal closure usually occurs to reduce the loss of water through transpiration, but at the same time, decreasing the entrance of CO 2 into leaves, limiting photosynthesis, and ultimately plant growth [3]. With an increase in the severity of drought, the functioning of photosynthesis can be further impaired by photochemical and biochemical dysfunctions occurring that was further crossed with C. arabica cv. Mundo Novo by IAC. The two genotypes display a relevant response ability to drought, although having different degrees of resilience as CL153 suffers a higher negative impact than Icatu in the photochemical and biochemical components of C-assimilation under severe drought [12,25]. Here, we explore the molecular mechanisms beyond such striking physiological and biochemical differences. Understanding the molecular mechanisms that ultimately determine the response of Coffea to climatic changes is crucial to mitigate their harmful effects and establish better scenarios to maintain the sustainability of the value chain of coffee.

Overview of the RNA-Seq Data from the Two Coffee Genotypes
RNA sequencing resulted in an average of 21.8 million paired-end reads per sample, generating an average of 16.9 and 19.3 million high-quality unique reads in Icatu and CL153, respectively (Table S2). A high proportion of reads was mapped to the corresponding reference genome since only an average of 16% and 14% of reads from Icatu and CL153, respectively, were not mapped. Statistical details for each replicate are depicted in Table S2.
The number of expressed genes varied in Icatu from 29199 (SWD-aCO 2 ) to 31041 (WW-eCO 2 ), while much lower values were found in CL153, ranging from 19558 (SWD-eCO 2 ) to 20463 (WW-aCO 2 ) ( Figure 1). A principal component analysis (PCA) based on gene expression generally clustered the samples according to the different water treatments ( Figure S1). In addition, MWD plants were usually closer to WW plants under eCO 2 but with SWD plants under aCO 2 conditions, especially when considering Icatu plants ( Figure S1). was further crossed with C. arabica cv. Mundo Novo by IAC. The two genotypes di a relevant response ability to drought, although having different degrees of resilien CL153 suffers a higher negative impact than Icatu in the photochemical and bioche components of C-assimilation under severe drought [12,25]. Here, we explore the m ular mechanisms beyond such striking physiological and biochemical differences. U standing the molecular mechanisms that ultimately determine the response of Cof climatic changes is crucial to mitigate their harmful effects and establish better scen to maintain the sustainability of the value chain of coffee.

Overview of the RNA-Seq Data from the Two Coffee Genotypes
RNA sequencing resulted in an average of 21.8 million paired-end reads per sam generating an average of 16.9 and 19.3 million high-quality unique reads in Icatu CL153, respectively (Table S2). A high proportion of reads was mapped to the corresp ing reference genome since only an average of 16% and 14% of reads from Icatu CL153, respectively, were not mapped. Statistical details for each replicate are depict Table S2.
The number of expressed genes varied in Icatu from 29199 (SWD-aCO2) to 3 (WW-eCO2), while much lower values were found in CL153, ranging from 19558 (S eCO2) to 20463 (WW-aCO2) ( Figure 1). A principal component analysis (PCA) base gene expression generally clustered the samples according to the different water ments ( Figure S1). In addition, MWD plants were usually closer to WW plants under but with SWD plants under aCO2 conditions, especially when considering Icatu p ( Figure S1).

Response of Differentially Expressed Genes (DEGs) to Drought and eCO2
The impact of MWD on the total number of DEGs was minimal under eCO2 in parison with plants under aCO2 (Figure 2A,B). The number of DEGs commonly trigg

Response of Differentially Expressed Genes (DEGs) to Drought and eCO 2
The impact of MWD on the total number of DEGs was minimal under eCO 2 in comparison with plants under aCO 2 (Figure 2A,B). The number of DEGs commonly triggered by both water deficits was much lower under eCO 2 than under aCO 2 . In fact, under MWD, a high number of DEGs were up-regulated under eCO 2 plants when compared to their aCO 2 counterparts, especially in Icatu. by both water deficits was much lower under eCO2 than under aCO2. In fact, under M a high number of DEGs were up-regulated under eCO2 plants when compared to aCO2 counterparts, especially in Icatu.  Table S3). The strong impact o SWD was also depicted on heatmaps visualizing the expression profile of all signif DEGs, where a high degree of variation was found under this harsher drought lev contrast with MWD, which promoted only minor effects ( Figure S2). The full list of D can be found in Tables S4 and S5. The majority of the up-and down-regulated DEGs found had no annotated func ( Figure 2C,D). The remaining ones were mostly involved in 'catalytic activities,' follo by 'binding' in the two genotypes regardless of [CO2] levels ( Figure 2C,D). Under a MWD triggered the down-regulation of DEGs involved in 'structural molecule act in Icatu. Instead, in CL153 'transporter activity' DEGs were down-regulated under but up-regulated under eCO2. Still, under MWD and eCO2 there was a down-regul of DEGs involved in 'transporter activity' in Icatu and 'transcription regulator activi CL153.
Under harsher drought conditions (SWD), Icatu aCO2-plants showed an up-re  Table S3). The strong impact of the SWD was also depicted on heatmaps visualizing the expression profile of all significant DEGs, where a high degree of variation was found under this harsher drought level in contrast with MWD, which promoted only minor effects ( Figure S2). The full list of DEGs can be found in Tables S4 and S5. The majority of the up-and down-regulated DEGs found had no annotated functions ( Figure 2C,D). The remaining ones were mostly involved in 'catalytic activities', followed by 'binding' in the two genotypes regardless of [CO 2 ] levels ( Figure 2C,D). Under aCO 2 , MWD triggered the down-regulation of DEGs involved in 'structural molecule activity' in Icatu. Instead, in CL153 'transporter activity' DEGs were down-regulated under aCO 2 but up-regulated under eCO 2 . Still, under MWD and eCO 2 there was a down-regulation of DEGs involved in 'transporter activity' in Icatu and 'transcription regulator activity' in CL153.
Under harsher drought conditions (SWD), Icatu aCO 2 -plants showed an up-regulation of DEGs involved in 'transporter activity', 'catalytic activity', and 'binding', but the latter two categories were involved in down-regulated DEG, together with molecular function regulators ( Figure 2C,D). CL153 plants under SWD and aCO 2 showed a high number of upand down-regulated DEGs associated with catalytic activity and binding functions. Under eCO 2 , while Icatu plants showed a down-regulation of DEGs involved in catalytic activities and binding, CL153 plants rather showed up-and down-regulated DEGs involved in these two categories, as well as in transcription regulator activities.

Drought and eCO 2 Impact on DEGs Associated with Specific Biochemical Pathways
Drought had a relevant impact on a high number of DEGs involved in respiration, antioxidant, and lipid biochemical pathways, especially in Icatu where a total of 342 DEGs (respiration: 168, antioxidant: 106, lipid metabolism: 68 DEGs; Table S6) were found to be affected, in comparison with 91 in CL153 (respiration: 31, antioxidant: 29, lipid metabolism: 31 DEGs; Table S7). Under MWD and aCO 2 , DEGs associated with these biochemical pathways were mostly down-regulated but they were substantially reduced under eCO 2 ( Figure 3). The positive effect of eCO 2 was even more relevant under SWD in Icatu plants, which showed an increase in up-regulated DEGs.  Table S7). Under MWD and aCO2, DEGs associated with these biochemical pathways were mostly down-regulated but they were substantially reduced under eCO2 ( Figure 3). The positive effect of eCO2 was even more relevant under SWD in Icatu plants, which showed an increase in up-regulated DEGs.  (Table S8) or CL153 (Table S9). The expression of the remaining photosynthetic-related DEGs decreased with drought with a few exceptions: transkelotase and RuBisCO activase 1 were up-regulated in Icatu plants in both water deficits (Table S8); in CL153 plants, the translocase was up-regulated in both water deficits but only under aCO2; while RuBisCO activase 1 was up-regulated under SWD and aCO2 but down-regulated under SWD and eCO2 (Table S9).  (Table S8) or CL153 (Table S9). The expression of the remaining photosynthetic-related DEGs decreased with drought with a few exceptions: transkelotase and RuBisCO activase 1 were up-regulated in Icatu plants in both water deficits (Table S8); in CL153 plants, the translocase was up-regulated in both water deficits but only under aCO 2 ; while RuBisCO activase 1 was up-regulated under SWD and aCO 2 but down-regulated under SWD and eCO 2 (Table S9).

DEGs Involved in the Response to Water Deprivation and Desiccation
Drought triggered a high number of DEGs involved in water deprivation and desiccation responses, especially in Icatu plants (155 DEGs) in comparison with CL153 (16 DEGs). The full list of DEGs is depicted in Table S10. DEGs involved in water-deprivation responses showed minor changes under MWD ( Figure 5). For instance, no significant DEGs were detected in CL153 plants grown under MWD and eCO 2 . Under SWD, most droughtresponsive genes were down-regulated under aCO 2 but up-regulated under eCO 2 , an effect mostly seen in Icatu plants ( Figure 5).
FOR PEER REVIEW 7 of 24

DEGs Involved in the Response to Water Deprivation and Desiccation
Drought triggered a high number of DEGs involved in water deprivation and desiccation responses, especially in Icatu plants (155 DEGs) in comparison with CL153 (16 DEGs). The full list of DEGs is depicted in Table S10. DEGs involved in water-deprivation responses showed minor changes under MWD ( Figure 5). For instance, no significant DEGs were detected in CL153 plants grown under MWD and eCO2. Under SWD, most drought-responsive genes were down-regulated under aCO2 but up-regulated under eCO2, an effect mostly seen in Icatu plants ( Figure 5).   Table 1 and Table S10). SWD and aCO 2 triggered the most down-regulation of the probable xyloglucan endotransglucosylase/hydrolase protein 6 (FC: −12.31) and the most up-regulation of the protein phosphatase 2C 51-like (FC: 5.62). Under eCO 2 and SWD, the lowest down-regulation was also found for the probable xyloglucan endotransglucosylase/hydrolase protein 6 (FC: −11.17), while the highest up-regulation was recorded for the Late Embryogenesis Abundant protein Dc3-like (LEA-DC3; FC: 7.19). Notably, a high number of DEGs were involved in the antioxidant system as the galactinol synthase 2-like, sucrose synthase 2-like, the homeobox-leucine zipper ATBH-12, as well as aquaporins were up-regulated in Icatu under SWD and eCO 2 (Table S10). By contrast, the minor responses recorded in CL153 plants and mostly recorded under SWD and aCO 2 involved a high number of Aspartic Protease in Guard Cell 1-like (ASPG1; 6 out of 16 total DEGs; Table 2). In this genotype, the 18 kDa seed maturation was the most up-regulated under MWD, independently of CO 2 levels (aCO 2 FC: 9.18; eCO 2 FC: 8.79) while the lowest was the APG1, also independently of CO 2 levels (aCO 2 FC: −7.93; eCO 2 FC: −8.16; Table 2). Remarkably, under aCO 2 , the dehydrin DH1a was also highly up-regulated under both water deficit levels (FC: 6.05 for MWD; FC 5.47 for SWD) while it was less expressed under eCO 2 and only under SWD (FC 2.37). In fact, under eCO 2 and MWD, no water-deprivation-related DEG was recorded. Under SWD, only one downregulated DEG was found (the putative movement binding protein 2C; FC:−19.92) while the remaining were up-regulated, especially the APG1 (FC: 3.96). Icatu plants under aCO 2 showed a moderate increase in enriched GO terms as water deficit severity increased, being all associated with down-regulated DEGs ( Figure 6).  Table S11). However, under SWD and eCO 2 , two enriched categories were linked to up-regulated DEGs: 'sequence-specific DNA binding' (GO:0043565) and the 'UDP−glycosyltransferase activity' (GO:0008194).
By contrast, drought triggered a higher number of enriched GO terms in CL153 than in Icatu (84 vs. 43) ( Figure 6; Table S12). Under aCO 2   Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG), only pathways linked to down-regulated DEGs were found to be significantly affected by drought (Figure 7). Under MWD and aCO 2 , the 'flavonoid biosynthesis' was the only KEGG pathway affected by drought and only in Icatu plants. Under SWD, the 'photosynthesis' pathway was mostly affected in Icatu and CL153 under aCO 2 , but only under eCO 2 in CL153 plants. In addition, the 'indole alkaloid biosynthesis' pathway was affected under eCO 2 in Icatu plants, while in CL153 plants, the 'photosynthesis antenna proteins' pathway was also affected under aCO 2 (Figure 7).
Based on the Kyoto Encyclopedia of Genes and Genomes (KEGG), only pathways linked to down-regulated DEGs were found to be significantly affected by drought (Figure 7). Under MWD and aCO2, the 'flavonoid biosynthesis' was the only KEGG pathway affected by drought and only in Icatu plants. Under SWD, the 'photosynthesis' pathway was mostly affected in Icatu and CL153 under aCO2, but only under eCO2 in CL153 plants.
In addition, the 'indole alkaloid biosynthesis' pathway was affected under eCO2 in Icatu plants, while in CL153 plants, the 'photosynthesis antenna proteins' pathway was also affected under aCO2 (Figure 7).

Differential Transcriptional Drought Regulation Responses in the Two Coffea Genotypes
Despite the information available, the molecular basis of stress responses in woody species is still limited in comparison with model plants. Understanding which genes may play a role under different degrees of drought severity, and which physiological traits are controlled by such genes, is of crucial importance for understanding the capabilities of trees to successfully cope with stress, especially in commercially important crops such as coffee. To maintain the global supply of this crop, it is important to promote the screening and development of tolerant varieties to face the increasingly expected impacts of drought events. In this study, we showed that the two genotypes, belonging to two different Coffea species, activated distinct transcriptional responses to cope with drought.
In comparison with CL153, Icatu showed (i) a higher number of expressed genes under MWD, but especially under SWD ( Figure 1); (ii) a higher number of specific DEGs, especially under SWD ( Figure 2); (iii) a higher number of DEGs involved in respiration, antioxidant activities, and lipid metabolism ( Figure 3); (iv) a higher activation of DEGs involved in light reactions of photosynthesis, the Calvin cycle and photorespiration (Figure 4); and (v) a higher number of DEGs involved in responses to water deprivation and desiccation ( Figure 5). This transcriptional result explains previous physiological studies that reported a higher photochemical performance in Icatu than in CL153, due to the preservation or reinforcement of photosynthetic components, strengthened enzymatic antioxidative system, and a greater abundance of photoprotective pigments [7,12,26,27].
In response to drought, CL153 rather activated a high number of aspartic proteases (ASPG1; Table 2). Aspartic proteases play a fundamental role in the response of plants to drought, especially during ROS increments [34]. ROS can act as a positive regulator of

Differential Transcriptional Drought Regulation Responses in the Two Coffea Genotypes
Despite the information available, the molecular basis of stress responses in woody species is still limited in comparison with model plants. Understanding which genes may play a role under different degrees of drought severity, and which physiological traits are controlled by such genes, is of crucial importance for understanding the capabilities of trees to successfully cope with stress, especially in commercially important crops such as coffee. To maintain the global supply of this crop, it is important to promote the screening and development of tolerant varieties to face the increasingly expected impacts of drought events. In this study, we showed that the two genotypes, belonging to two different Coffea species, activated distinct transcriptional responses to cope with drought.
In comparison with CL153, Icatu showed (i) a higher number of expressed genes under MWD, but especially under SWD ( Figure 1); (ii) a higher number of specific DEGs, especially under SWD ( Figure 2); (iii) a higher number of DEGs involved in respiration, antioxidant activities, and lipid metabolism ( Figure 3); (iv) a higher activation of DEGs involved in light reactions of photosynthesis, the Calvin cycle and photorespiration ( Figure 4); and (v) a higher number of DEGs involved in responses to water deprivation and desiccation ( Figure 5). This transcriptional result explains previous physiological studies that reported a higher photochemical performance in Icatu than in CL153, due to the preservation or reinforcement of photosynthetic components, strengthened enzymatic antioxidative system, and a greater abundance of photoprotective pigments [7,12,26,27].
In response to drought, CL153 rather activated a high number of aspartic proteases (ASPG1; Table 2). Aspartic proteases play a fundamental role in the response of plants to drought, especially during ROS increments [34]. ROS can act as a positive regulator of ABA signaling in guard cells, but the excessive accumulation of ROS during stress can be toxic [35]. Thus, to balance ROS production and scavenging, either the levels are modulated through signal transduction, or the cells have mechanisms to detoxify excessive ROS during stress [36]. Hence, the up-regulation of aspartic proteases in CL153 may help to scavenge the excessive amount of ROS. The up-regulation of dehydrins (namely DH1a) in both water deficits may also be involved in protective reactions to dehydration [37]. Dehydrins were also found to be highly expressed in the leaves of drought-stressed Coffea plants such as C. arabica cvs. Catuaí and Mundo Novo, C. canephora cv. Apoatã [38], as well as in the genotypes studied here [39], where they seem to play a key role in the acclimation response of Coffea.
Our transcriptomic results also found a high number of DEGs involved in the antioxidant system, especially in Icatu (Figure 3). This suggests the action of three mechanisms involved in drought acclimation: (i) activation of antioxidant activities to scavenge the excessive amount of ROS and reduce oxidative damage in plants, including enzymes (e.g., catalase, superoxide dismutases, peroxidases) and non-enzymatic molecules (e.g., ascorbic acid, α-tocopherol, raffinose family oligosaccharides), in agreement with studies from other plants [40,41]; (ii) antioxidative mechanisms complemented with thermal dissipation mechanisms (e.g., photoprotective carotenoids) and changes in the cyclic electron flow (CEF) to protect the photosystem (PS) I and/or II. In fact, (i) and (ii) are transversally found in resilient Coffea genotypes as response mechanisms to drought [12,26], heat [42][43][44], cold [45,46], and high irradiance [47,48] stresses. Additionally, the activation of other molecules, such as aquaporins reported here in Icatu and in previous studies involving other C. arabica genotypes [49] or the heat shock protein 70 kDa as reported also in other Coffea studies [39,50]. Altogether, results reveal the expression of DEGs, whose products are important to limit or control water loss in Coffea, regulating stomatal closure in leaves subjected to drought conditions. Their overexpression upon stress was validated by qRT-PCRs (Figure 8), supporting the action of several antioxidant molecules in response to drought. It also corroborates the results of RNA-sequencing, providing a set of candidate genes involved in drought tolerance in coffee plants. These mainly include genes encoding antioxidant enzymes involved in the biosynthesis of small antioxidant molecules and water and ion movement, such as aquaporins and ion transporters or the biosynthesis of osmolytes. Other proteins that function directly in the protection of proteins and membranes, such as late embryogenesis abundant proteins, heat shock proteins, protein phosphatases, and aspartic proteases also play important roles in Coffea resilience to abiotic stresses. This vast information will help in the identification of target sites suitable for gene editing, which will undoubtfully make progress in the future.

Influence of eCO 2 in Drought-Responsive Genes
The level and type of transcripts found in this study support the positive effects of eCO 2 in mitigating some impacts of drought, especially in the functioning and components of the photosynthetic apparatus [12]. The up-regulation of antioxidant-related DEGs under eCO 2 , mostly under SWD and in Icatu plants (Figure 3) is congruent with physiological measurements that showed an increase in Cu,Zn-SOD activity under MWD, rising further under SWD, but only under eCO 2 [27]. The up-regulation of antioxidant-related DEGs under eCO 2 might strengthen the scavenging potential for molecular O 2 when the photochemical use of energy is small, reducing the probability of electron flow to lower energy states [12], helping to alleviate some impacts of drought in the photosynthetic machinery ( Figure 4).
The imposition of eCO 2 to SWD plants promoted the up-regulation of several important stress response genes, including the LEA-DC3 in Icatu and the APG1 in CL153 as mentioned before. Late embryogenesis abundant proteins compose the most abundant and characterized group of intrinsically disordered proteins that prevent and repair the damage caused by environmental stresses. A positive association between the accumulation of LEA and environmental stresses, such as drought, heat, and salinity, has been outlined in several other plant species, where they bind to enzymes to prevent the loss of activity under stressful conditions [51,52]. This is consistent with the biological functions of LEAs namely in oxidant scavenging activities, enzyme and nucleic acid preservations, and the membrane maintenance that occurs in genotypes/species that can better cope with environmental stresses [53].
Directly or indirectly, ABA signaling pathways, which include ABA-dependent, ABAresponsive element/ABRE-binding factors (ABRE/ABF), as well as ABA-independent genes as dehydration-responsive elements regulate Coffea responses to stress. ABA is involved in both drought-induced and eCO 2 -induced stomatal closure in a dual way, including root-derived and foliar ABA. For instance, under well-watered conditions, C. arabica plants grown under eCO 2 showed lower whole-plant transpiration rates than under aCO 2 [21]. These changes, although unrelated to stomatal conductance or foliar ABA levels, are associated with faster stomata closure rates upon rapid increases in vapor pressure deficit under eCO 2 [21]. For instance, during exposure to drought, Coffea plants grown under eCO 2 can maintain higher water potentials and plant hydraulic conductance than under aCO 2, due to a higher transcript abundance of aquaporins [21]. In the genotypes studied here, drought alone prompted gradual ABA increases of~46% in MWD plants in both genotypes, and 100% (CL 153) and 184% (Icatu) under SWD conditions, whereas single eCO 2 increased ABA levels (by 85%) but only in Icatu [12]. This is important since ABA controls stomatal closure under stress, which is one of the first defense responses to reduce water loss by restricting the transpiration flow in response to a rising air evaporative demand or a decreased soil water availability [54]. This suggests that, at least under MWD, eCO 2 seemed to decouple ABA action from stomatal closure in both genotypes. In agreement with this hypothesis, dehydration is postponed in MWD plants under eCO 2 , and a delayed stomata response to soil drying under eCO 2 is found in some coffee genotypes [12,21,22], maintaining greater stomatal conductance than expected based on ABA concentrations. Such a greater stomatal opening under MWD would allow greater C-assimilation gains under eCO 2 [12], also revealing a more profound impact of eCO 2 than the direct stimulation of C-assimilation. A better understanding of ABA concentrations in the xylem sap is necessary to understand the sensitivity of Coffea stomata to [ABA] xylem.

Overall Regulatory Mechanisms Involved in Coffea Responses
MWD had minor impacts under aCO 2 with few enriched categories being found in Icatu, and all were associated with down-regulated DEGs, while CL153 plants showed a high increase in GO categories, and were mostly associated with up-regulated DEGs ( Figure 6).
Under SWD, a few enriched categories were even up-regulated in Icatu ('sequencespecific DNA binding' and the 'UDP−glycosyltransferase activity') and CL153 ('defense response', 'regulation of transcription DNA-templated', 'oxidoreductase activity' and 'transcription regulator activity'). Genes assigned to these GO terms are usually involved in a high number of developmental processes and stress responses, plant hormone activation, and the production of antioxidants in response to stresses, including drought [55]. In Coffea, as in other species, glycosylation catalyzed by glycosyltransferases, as well as oxidoreductase activities, play an essential role in regulating the stability, availability, and biological activity of antioxidant compounds and the integrity of cellular membranes [56]. This is crucial to cope with the effects of water deprivation in Coffea since tolerance is associated with the ability of tissues to withstand low water potentials and plant membrane transport systems play a significant role under water scarcity [57]. Depending on energy needs, translocation through biological membranes occurs, passively or actively, but drought can compromise the integrity of membranes and is, therefore, relevant an up-regulation of categories linked to cellular membranes in the two Coffea genotypes, even under MWD.
The 'flavonoid biosynthesis' KEGG pathway was affected in Icatu plants grown under MWD and aCO 2 (Figure 7). Plant phenolic compounds, especially flavonoids, can provide resistance to biotic and abiotic stresses, and can be enhanced upon drought [58]. As nonenzymatic antioxidants, hydroxyl groups in flavonoids participate in the scavenging of oxygen free radicals, alleviating stress-induced oxidative damage as is widely reported [59], thus, affecting acclimation responses of Icatu, at least under aCO 2 . However, SWD affected the KEGG photosynthetic pathway, namely in Icatu plants under aCO 2 , as well as in CL153 plants under both [CO 2 ] levels together with the 'photosynthesis antenna proteins' pathways in CL153 plants under SWD and aCO 2 (Figure 7). Since these pathways are both linked to down-regulated DEGs, this would imply some effects in the photosynthetic machinery of these genotypes, especially in CL153. Indeed, a large dilution of the impacts of drought on the net photosynthesis in these genotypes was found to be promoted by eCO 2 under MWD, consistent with a tendency to the maintenance of the PSII efficiency and higher PSs activity in both genotypes while under SWD, the net photosynthesis and stomatal conductance were severely reduced, regardless of [CO 2 ] [12]. Nevertheless, even under SWD, a relevant potential for C-assimilation was preserved, with the photosynthetic capacity (A max ) showing values close to 60% (CL153), or even higher than 70% (Icatu) relative to those displayed by their respective WW controls [12]. This photochemical protective mechanism results in a lower need for dissipation processes and a reduced PSII inhibition status [43] consistent with the down-regulation of photosynthetic-related DEGs found in this study (Figure 4) or the KEGG enrichment results (Figure 7) since no new molecules would be needed due to the protective mechanisms involved.
A minor inhibition in the 'indole alkaloid biosynthesis' (IAB) KEGG pathway was also recorded in Icatu under SWD and eCO 2 . Plant peroxidases may accept alkaloids as substrates, as well as phenols and flavonoids [60], and metabolize H 2 O 2 as an electron donor for phenol peroxidases, resulting in the formation of phenoxyl radicals, which can be regenerated by a non-enzymatic reaction with an ascorbate function as an H 2 O 2 scavenging system [61]. Thus, either the IAB KEGG pathway is not involved/not needed in acclimation responses of Icatu or this would suggest a decrease in the defense ability of this genotype against oxidative stress. However, that contrasts with the high activity found in protective molecules and antioxidant enzymes in Icatu [12,27], and the large reinforcement of Cu,Zn-SOD, APX, and catalase (CAT) activities [7].

Evidence of Post-Transcriptional Regulatory Mechanisms in Coffea Responses to Stress
Previous findings showed that Coffea plants, namely Icatu, can maintain the potential photosynthetic functioning under the imposition of SWD due to a greater antioxidative response, which contrasts with the transcriptomic results shown here where photosyntheticrelated DEGs were mostly down-regulated (Figures 3 and 4) and the KEGG photosynthetic pathway was highly affected under SWD (Figure 7). Physiological studies showed that increasing drought severity progressively affected gas exchange and fluorescence parameters in both genotypes, with non-stomatal limitations becoming gradually dominating, and having strong impacts on the photochemical and biochemical components and functioning of Coffea, especially in CL153 plants under SWD and aCO 2 [12]. In contrast, Icatu plants were tolerant to SWD, with minor, if any, negative impacts on the potential photosynthetic functioning and components, e.g., A max , F v /F m , electron carriers, photosystems (PSs) and RuBisCO activity, under aCO 2 [12]. Under MWD, eCO 2 delayed stress severity and promoted photosynthetic functioning in both genotypes, with lower energy dissipation, while stomatal closure was decoupled from increases in ABA. Under SWD, most of the negative impacts felt on the photosynthetic components and their potential performance were reduced under eCO 2 , at least considering CL153, since Icatu was barely affected in both [CO 2 ] levels under SWD [12]. Still, strong effects were detected in RuBisCO, as the most sensitive photosynthetic component [12]. However, proteomic analyses have also shown a higher abundance of drought-responsive proteins in Icatu than in CL153, together with enriched GO terms, and enriched KEGG pathways associated with stress responses and the control of oxidative stress categories found here [12,27,62]. Thus, these contrasting results suggest the existence of important post-transcriptional regulation in Coffea, at least in the genotypes studied here. Other studies have also highlighted that protective metabolites often do not show a clear pattern between transcript accumulation and metabolite/physiological responses in response to stress, being likely to not be transcriptionally regulated. For instance, [63] studied the transcript response to eCO 2 in Solanum lycopersicum and its wild relative S. pennelli, and no clear transcriptomic pattern was found, but rather a translational regulatory mechanism, hypothetically involved in the differential ribosomal loading of transcripts in the two species. Additionally, in Saccharomyces cerevisiae, relatively few (~15%) of the mRNAs that were translationally up-regulated in response to H 2 O 2 showed similar increases in transcript levels [64], revealing a complex transcriptional and translational reprogramming to stress. The existence of a complex translational program in Coffea would also explain the physiological and biochemical performance of these genotypes [12,26], as well the amplified acclimation responses at the proteomic level [27], namely in Icatu plants, despite the down-regulation of transcripts reported here.

Imposition and Monitoring of Water Deficit Conditions and Sampling
Plants previously maintained without water restriction were divided into three groups. The first one was maintained under well-watered (WW) conditions, with a leaf predawn water potential (Ψ pd ) above −0.35 MPa. In the other two groups, drought was imposed by a gradual reduction of irrigation, allowing plants to express their potential acclimation ability for two weeks, to promote Ψ pd decline to values between −1.5 and −2.5 MPa (moderate water deficit-MWD) or below −3.5 MPa (severe water deficit-SWD), representing ca. 80 (WW), 35 (MWD) and 10% (SWD) of maximal water availability in pots [6]. These MWD and SWD conditions were maintained for another two weeks by adding adequate water amounts according to each water deficit level. Samples were then collected for transcriptomic analysis. Exceptionally, Icatu eCO 2 plants under MWD were submitted to total water withholding in the last 5 days of the 4-week period, to force the reduction of Ψ pd , which, even so, did not shift below −0.6 MPa. Leaf Ψ pd was determined immediately after leaf excision, using a pressure chamber (Model 1000, PMS Instrument Co., Albany, OR, USA).

RNA Extraction and Illumina Sequencing
Newly matured leaves from plagiotropic and orthotropic branches from the upper third part (well illuminated) of each plant were collected under photosynthetic steady-state conditions after 2 h of illumination, flash frozen in liquid nitrogen, and stored at −80 • C. Total RNA was extracted from 36 samples (two genotypes × three water treatments × two [CO 2 ] × three individual plants) using the Analytik-Jena InnuSPEED Plant RNA Kit (Analytik Jena Innuscreen GmbH, Jena, Germany) following [50]. RNA quantity and quality were determined using a BioDrop Cuvette (BioDrop, UK) and an Agilent 2100 Bioanalyzer

Quality Analysis of Sequencing Data
Raw reads were processed using FastQC version 0.11.9 [66] to remove low-quality reads. FastQ Screen version 0.14 [67] was used to check for contaminants against the genome of the most common model organisms (e.g., Homo sapiens, Mus musculus, Rattus norvegicus, Drosophila melanogaster, Caenorhabditis elegans, Saccharomyces cerevisiae, Escherichia coli) and adapter databases (e.g., Mitochondria RNA, PhiX, Vector from UniVec database, FastQ Screen rRNA custom database, and FastQ Screen Adapters database). Since all reads presented an overall good quality, the trimming step was skipped. Recent studies showed that this process is redundant in the quantification of expression data from RNA-seq since most aligners can perform soft-clipping to effectively remove adapter sequences and rescue low-sequencing-quality bases that would be removed by read trimming tools, improving the accuracy in the quantification of gene expression [68].

Reference-Based Mapping and Assembly
The raw reads of Icatu were mapped to the reference genome of C. arabica downloaded from the NCBI (https://www.ncbi.nlm.nih.gov/assembly/GCF_003713225.1, accessed on 4 April 2021), while the raw reads of CL153 were mapped to the C. canephora genome downloaded from the Coffee Genome Hub (http://coffee-genome.org/download, accessed on 4 April 2021) [69] using STAR version 2.7.8a [70]. Htseq-count v0.11.0 [71] was then used to quantify only uniquely mapped genes. Samtools version 1.12 [72] and gffread version 0.12.1 [73] were used throughout the analysis to obtain general statistics on genome mapping. A principal component analysis (PCA) was performed on the expression data of genes, TMM (trimmed mean of means) normalized and log10-transformed using ggplot2 version 3.3.3 library [74] of R software version 4.0.2 [75]. Through visual inspection of the PCA, replicate 7B was considered an outlier and excluded from downstream analyses ( Figure S1).

Identification of Differentially Expressed Genes (DEGs)
DEGs from the plants grown under different water treatments and different [CO 2 ] were estimated and compared as follows: MWD-aCO 2 vs. WWa-CO 2 , SWD-aCO 2 vs. WW-aCO 2 , MWD-eCO 2 vs. WW-eCO 2 , and SWD-eCO 2 vs. WW-eCO 2 . Differential expression analyses were performed using the DEGs commonly found by both DESeq2 version 3.8 [76] and edgeR version 3.26.0 [77]. The resulting values of expression were adjusted using Benjamini and Hochberg's approach to control the false discovery rate (FDR; [78]). Genes with a normalized non-zero log 2 fold change (FC) expression and an FDR < 0.01 in both tools were defined as differentially expressed. Python's matplotlib library was used to plot Venn diagrams and bar plots [79]. Using ggplot2, heatmaps with dendrograms were plotted to visualize DEGs based on the differential expression patterns between the different comparisons. To prevent high DEGs from clustering together without considering their expression pattern, log 2 FC was scaled by gene expression across treatments (row Z-score).

Drought and eCO 2 Impact on DEGs Associated with Specific Biochemical Pathways
Due to the fundamental role of photosynthesis, respiration, lipid profile changes, and the antioxidant system in the process of coffee acclimation to environmental stresses, a specific/fine-tuned search was performed among the significant DEGs associated with these processes. According to the reference genome and the UniProtKB database DEGs annotated with the following direct and child GO terms were searched to study their regulation pattern:

DEGs Involved in the Response to Water Deprivation and Desiccation
To better understand the impacts of water deficit, a specific search was performed, as described previously, among DEGs annotated with 'water transport' (GO:0006833), 'water homeostasis' (GO:0030104), 'response to water' (GO:0009415), 'response to water deprivation' (GO:0009414), and 'response to desiccation' (GO:0009269).

Functional Classification of Responsive DEGs
DEGs from CL153 and Icatu comparisons were annotated following the functional annotation of the reference genomes of C. canephora and C. arabica, respectively, as stated previously. GO enrichment analyses were applied to understand the functional classification of responsive DEGs through another over-representation analysis (ORA) using gProfiler [81] under g:SCS < 0.01. Results were summarized using REVIGO [82] by removing redundant GO terms within a similarity = 0.5. Enrichment non-redundant results were plotted using ggplot2 version 3.3.3 library, using the number of DEGs annotated with each term to set a Counts > 10 cut-off [74]. Since Coffea genomic annotations are not complete, namely in terms of KEGG pathways, DEGs were mapped to their Arabidopsis thaliana homologs against a local Swissprot database, filtering gene hits by maximum e-value of 1.0 × 10 −3 and minimum identities of 40% [83], and using blastx from the Basic Local Alignment Search Tool (BLAST) version 2.10.1 command-line tool from the NCBI C++ Toolkit. These annotations were then used to perform an over-representation analysis (ORA) with gProfiler, searching for significantly (g:SCS < 0.01) enriched KEGG pathways.

Quantitative RT-PCR
Twelve transcripts were randomly selected for real-time quantitative PCR (qRT-PCR) to verify the accuracy of the levels of expression obtained under RNA-seq. Genes included: GMPM1: 18 kDa seed maturation; PP2C-51: protein phosphatase 2C 51-like; LEA-DC3: late embryogenesis abundant protein Dc3-like; DH1a: dehydrin DH1a; ATHB22: homeobox leucine zipper; SUS2: sucrose synthase 2-like; PIP2-2: aquaporin PIP2-2-like; XTH6: xyloglucan endotransglucosylase/hydrolase protein 6; GOLS2: galactinol synthase 2-like; CuSOD1: Superoxide dismutase [Cu-Zn]; APX Chl : chloroplast ascorbate peroxidase. All primer sequences are presented in Table S1. The primers were designed using Primer3 web version 4.1.0 [84] with an e-value < 2 × 10 −4 and a score >41. cDNA was synthesized from 1 µg total RNA using the SensiFASTTM cDNA Synthesis kit (Meridian BioScience, Cincinnati, OH, USA), according to the manufacturer's recommendations. The presence of a single amplification product of the expected gene size was verified by electrophoresis on a 1.5% agarose gel. PCR reactions were prepared using the SensiFASTTM SYBR No-ROX kit (Meridian BioScience, USA) according to the manufacturer's protocol. One negative sample was included for each primer pair, in which cDNA was replaced by water. Reactions were carried out in 96-well plates using a qTOWER 2.2 Thermal Cycler (Analytik, Jena, Germany) using the following parameters: hot start activation of the Taq DNA polymerase at 95 • C for 10 min, followed by 40 cycles of denaturation at 95 • C for 15 s, annealing at 60 • C for 30 s, elongation at 72 • C for 30 s. A melting curve analysis was performed at the end of the PCR run by a continuous fluorescence measurement from 55 • C to 95 • C with sequential steps of 0.5 • C for 15 s. A single peak was obtained and no signal was detected in the negative controls. Three technical replicates were used for each analyzed plant. Gene expression was quantified using malate dehydrogenase (MDH) and ubiquitin (UBQ10) as reference genes [85]. To understand the agreement of these results with the one from RNA-sequencing, heatmaps were constructed considering the levels of transcripts and their expression levels from qRT-PCRs after being log 2 FC scaled by gene expression across treatments.

Conclusions
In this study, we showed how the single and combined effects of drought and eCO 2 triggered wide but differential responses at the transcriptional level in Coffea.
MWD had a minor impact on the number of transcripts differentially regulated by Icatu and CL153, contrary to SWD where a high number of DEGs were reported, being mostly down-regulated. eCO 2 attenuated the impacts of drought in the two genotypes, but especially in Icatu, in agreement with the contrasting physiological tolerance previously reported in these genotypes.
There was a predominance of protective and ROS-scavenging genes, directly or indirectly related to ABA signaling pathways involving Coffea tolerance responses. These genes were also involved in water deprivation and desiccation processes, such as LEA and protein phosphatases in Icatu and Aspartic Protease in Guard Cell 1-like and dehydrins in CL153, being their expression confirmed by qRT-PCR.
Enrichment analysis of GO and KEGG pathways revealed different regulatory mechanisms of Icatu and CL153 in response to drought, agreeing with the minor effects of MWD and the positive action of eCO 2 . However, a clear effect on photosynthetic pathways was recorded, namely under SWD and eCO 2 , contrary to previous physiological and biochemical studies.
The existence of a complex post-transcriptional regulatory mechanism is suggested to occur in Coffea explaining the discrepancies between transcriptional vs. proteomic and physiological data in these genotypes.