Understanding the Impact of Drought in Coffea Genotypes: Transcriptomic Analysis Supports a Common High Resilience to Moderate Water Deﬁcit but a Genotype Dependent Sensitivity to Severe Water Deﬁcit

: Water scarcity is the most signiﬁcant factor limiting coffee production, although some cultivars can still have important drought tolerance. This study analyzed leaf transcriptomes of two coffee cultivars with contrasting physiological responses, Coffea canephora cv. CL153 and Coffea. arabica cv. Icatu, subjected to moderate (MWD) or severe water deﬁcits (SWD). We found that MWD had a low impact compared with SWD, where 10% of all genes in Icatu and 17% in CL153 reacted to drought, being mainly down-regulated upon stress. Drought triggered a genotype-speciﬁc response involving the up-regulation of reticuline oxidase genes in CL153 and heat shock proteins in Icatu. Responsiveness to drought also included desiccation protectant genes, but primarily, aspartic proteases, especially in CL153. A total of 83 Transcription Factors were found engaged in response to drought, mainly up-regulated, especially under SWD. Together with the enrollment of 49 phosphatases and 272 protein kinases, results suggest the involvement of ABA-signaling processes in drought acclimation. The integration of these ﬁndings with complementing physiological and biochemical studies reveals that both genotypes are more resilient to moderate drought than previously thought and suggests the existence of post-transcriptional mechanisms modulating the response to drought.


Introduction
Along with the rapid expansion in population and global warming, water scarcity has become a worldwide challenge for agriculture [1][2][3]. To cope with water deficits, plants trigger a wide range of responses at the molecular, biochemical, and physiological

Water Stress Imposition and Leaf Water Status
Water conditions were imposed as previously described [25]. Briefly, the plants were divided into three groups. In the first group, individuals were maintained well irrigated (WW) throughout the experiment. In the other two groups, water deficit was gradually imposed during two weeks by partially withholding irrigation (with a partial water replacement of the amount lost in each pot) until stability of predawn leaf water potential (Ψpd) to plant values between −1.5 and −2.5 MPa (moderate water deficit-MWD) or below −3.5 MPa (severe water deficit-SWD). WW plants were maintained under full irrigation (Ψpd > −0.35 MPa). These conditions represented approximately 80% (WW), 25% (MWD) or 10% (SWD) of maximal water availability in pots [27]. After reaching the desired Ψpd values for MWD or SWD conditions, the pot moisture was maintained for another two weeks before leaf sampling. Leaf Ψpd was determined at predawn 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 steadystate conditions after 2 h of illumination, flash-frozen in liquid nitrogen and stored at −80 • C, being finely powdered in liquid N 2 prior to analysis. Total RNA was extracted from 18 samples (two genotypes × three water treatments × three biological replicates) using 20 mg of the frozen leaves. RNA was extracted using the Analytik-Jena InnuSPEED Plant RNA Kit (Analytik Jena Innuscreen GmbH, Berlin, Germany) following [34]. RNA quantity and quality were determined using BioDrop Cuvette (BioDrop, Manchester, UK) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). RNA Integrity Number (RIN) for the samples ranged from 8.89 to 9.22. The messenger RNA (mRNA) libraries were constructed with the Illumina "TruSeq Stranded mRNA Sample Preparation kit" (Illumina, San Diego, CA, USA) and sequenced on an Illumina NovaSeq6000 at Macrogen facilities (Macrogen, Geumcheon-gu, Seoul, Korea). Raw reads have been deposited in the NCBI Sequence Read Archive, BioProject accession PRJNA729673.

Quality Analysis of Sequencing Data
High-quality reads were obtained after several steps of quality checks, including trimming and removal of adaptor/primer and low-quality reads using FastQC version 0.11.8 [35] and Trimmomatic version 0.38 [36]. FastQ Screen version 0.13 [37] was used to check for contaminants against the genome of the most common model organisms and adapter databases (e.g., Mitochondria RNA, PhiX, Vector from UniVec database, FastQ Screen rRNA custom database and FastQ Screen Adapters database).

Reference-Based Mapping and Assembly
The filtered high-quality reads were mapped to the reference genome of C. canephora downloaded from the Coffee Genome Hub (http://coffee-genome.org/download, accessed on 4 January 2020) [38] using STAR version 2.6.1 [39]. Htseq-count v0.11.0 [40] was used to count uniquely mapped genes. Samtools version 1.9 [41] and gffread version 0.9.9 [42] were used throughout the analysis to obtain general statistics of the genome mapping.

Identification of Differentially Expressed Genes (DEGs)
Gene expression normalization of all samples was estimated with the DESeq method (median of ratios normalization) to account for sequencing depth and RNA composition, which is appropriate for the differential expression (DE) analysis across samples. A Principal Component Analysis (PCA) was performed on the expression data of genes, FKPM normalized and log10-transformed, using the ggplot2 version 3.3.2 library [43] of R software version 3.5.1 [44]. Through visual inspection of the PCA plot, the 7B replicate was considered an outlier and thus excluded from the downstream analysis ( Figure S1). DESeq2 v1.28.1 [45] was used to check for differences in the relative abundance of the genes between the different water treatments. The Benjamini-Hochberg approach was used for controlling the false discovery rate, FDR [46]. Differentially expressed genes (DEGs) were defined as genes with a normalized non-zero log2 fold change (FC) expression and an FDR < 0.01. Python's matplotlib library was used to plot Venn diagrams and bar plots [47]. The functional annotation of the reference genome of C. canephora referred above was used to search the top responsive DEGs. To better understand the effects of drought, a specific analysis was performed among the DEGs annotated with the direct and child GO terms "response to water deprivation" (GO:0009414) and "response to desiccation" (GO:0009269). Additionally, due to its crucial importance in the acclimation response of coffee plants to changing environmental conditions [23,33,48], a specific search was performed among DEGs annotated with the following terms: "antioxidant activity" (GO:0016209), "response to oxidative stress" (GO:0006979) under Antioxidant activity; "cellular respiration" (GO:0045333), "mitochondrion" (GO:0005739), "malate dehydrogenase activity" (GO:0016615), "pyruvate kinase activity" (GO:0004743) under Cellular respiration; "fatty acid metabolic process" (GO:0006631) and LOX (GO:0004051, GO:0016702) under Lipid metabolism; and "photosynthesis" (GO:0015979), "photosystem" (GO:0009521), "photosynthetic electron transport chain" (GO:0009767), "photorespiration" (GO:0009853) and "chlorophyll biosynthetic process" (GO:0015995) under Photosynthesis.

Regulation Patterns of Transcription Factors
Blastx from the Basic Local Alignment Search Tool (BLAST) version 2.10.1 command line tools from the NCBI C++ Toolkit was used to map the DEGs against Arabidopsis thaliana homologs using a local Swissprot database, filtering gene hits by maximum E-Value of 1.0 × 10 −3 and minimum identities of 40% [49]. Then, to study the regulation pattern of transcription factors (TFs) among the detected DEGs, a list of A. thaliana TFs related to drought was retrieved from DroughtDB (http://pgsb.helmholtz-muenchen. de/droughtdb/, accessed on 15 May 2021) and searched among DEGs. To complement these analyses, TFs were also searched among DEGs if annotated with "DNA-binding transcription factor activity" (GO:0003700) and "general transcription initiation factor binding" (GO:0140296) in the reference genome. To understand the enrollment of protein kinases and phosphatases in the regulation of drought, the following GO terms were also searched among DEGs: "phosphatase activity" (GO:0016791) and "protein kinase activity" (GO:0004672).

Enrichment Analysis of Gene Ontology
Gene Ontology (GO) enrichment analyses were applied to understand the functional classification of DEGs through an over-representation analysis (ORA), using gProfiler [50] under FDR < 0.05. REVIGO [51] was used to summarize results by removing redundant GO terms with a similarity ≥0.5. Enrichment nonredundant data with these FDR and similarity cutoffs were plotted using ggplot2. This same package was used to plot heatmap with dendrograms to visualize DEGs based on the differential expression patterns between the different treatments. To prevent highly differentially expressed genes from clustering together without considering their expression pattern, log2 fold change was scaled by gene across treatments (row Z-score).

Overall Transcriptome Profiling and Mapping Statistics
Quality analysis, data trimming, and filtering generated an average of 24.8 million (CL153) and 22.7 million (Icatu) clean reads, from an average of 25.1 and 23.0 million raw reads, respectively. Overall, an average of 90% and 84% cleaned reads from CL153 and Icatu, respectively, were mapped to the reference genome (Table S1). Statistical details for each replicate are depicted in Table S1.
Besides the specificity linked with the level of water deficit, results also showed that the two genotypes reacted differently to drought as only a small number of DEGs were commonly found between genotypes and between water deficit treatments (115 downregulated and 27 up-regulated) ( Figure 2C). Both CL153 and Icatu showed a low number  Figure 2C). This reveals a high degree of up-regulated DEGs in CL153 under SWD, whereas, in Icatu, the number of specific down-regulated DEGs was 2.4 times higher than the up-regulated ones. The same pattern was found in the heatmap of treatment-specific DEGs that showed a small genotype differentiation under MWD (Figure 3). In sharp contrast, a higher degree of variation was found under SWD, especially in CL153, where DEGs were more up-regulated in this genotype than in Icatu, in agreement with the previous analysis. Besides the specificity linked with the level of water deficit, results also showed tha the two genotypes reacted differently to drought as only a small number of DEGs wer commonly found between genotypes and between water deficit treatments (115 down regulated and 27 up-regulated) ( Figure 2C Down-regulated Up-regulated 57 63 A.

B.
C. genotype differentiation under MWD ( Figure 3). In sharp contrast, a higher degree of variation was found under SWD, especially in CL153, where DEGs were more upregulated in this genotype than in Icatu, in agreement with the previous analysis.

Identification and Classification of DEGs
Under drought, the 10 top up-regulated DEGs in CL153 were primarily involved in oxidoreductase activity and FAD binding. For instance, under MWD, top up-regulated DEGs in CL153 included several reticuline oxidase-like genes, as well as different cytochrome P450 genes, a carotenoid cleavage dioxygenase, and the TF ERF027 (Tables 1  and S3).

Identification and Classification of DEGs
Under drought, the 10 top up-regulated DEGs in CL153 were primarily involved in oxidoreductase activity and FAD binding. For instance, under MWD, top up-regulated DEGs in CL153 included several reticuline oxidase-like genes, as well as different cytochrome P450 genes, a carotenoid cleavage dioxygenase, and the TF ERF027 (Tables 1 and S3).
Under SWD, the top DEGs in CL153 were involved in similar functions and at the same level of fold changes as under MWD, also showing an up-regulation of reticuline oxidase genes and mostly of an acid phosphatase gene (PAP20) involved in hydrolase activity and metal ion binding (Tables 2 and S4). However, while under MWD, CL153 showed a downregulation of DEGs mostly involved in binding, auxin production, and transporter activity, under SWD, the effect was two times higher, showing a strong down-regulation of the ROP-interactive CRIB motif-containing protein 4 gene (FC = −21.55; Table 2). Under MWD, several heat shock proteins were among the top 10 up-regulated DEGs in Icatu, namely the glycoside hydrolase family 79 gene that showed the highest regulation (FC = 21.23), while the TF ORG2 was the most down-regulated DEG (FC = −25.66; Tables 3 and S5).
Under SWD, Icatu top up-regulated DEGs were involved in binding and transporter activity but mainly on transferase activities involving the UDP-glycotransferase 75D1 (FC = 20.01) while down-regulating the TF ORG2 as reported under MWD (Tables 4 and S6).

Regulation Patterns of DEGs Directly Linked to Water Deprivation and Desiccation
To better understand the impacts of water deficit, a specific search performed among DEGs annotated with "response to water deprivation" (GO:0009414) or "response to desiccation" (GO:0009269) found 22 additional DEGs, mostly expressed under SWD in the two genotypes (Table 5). In CL153, these drought-responsive DEGs were slightly upregulated in MWD (5 out of 8) and more down-regulated under SWD (8 out of 14). In Icatu, these DEGs were all down-regulated under MWD (0 up and 7 down) and mostly under SWD (3 up and 9 down). Notably, the Desiccation protectant protein Lea14 was found to be commonly up-regulated by CL153 (MWD and SWD) and Icatu (only under SWD). Additionally, a large majority of DEGs (10 out of 22) were linked to the Aspartic Protease in Guard Cell 1 gene (ASPG1), being all down-regulated in Icatu under the two water deficits, while some were up-regulated in CL153. Many of these drought-responsiveness DEGs were linked with the 4th chromosome of C. canephora ( Figure S3). Table 5. Differentially expressed genes (DEGs) under moderate water deficit (MWD) or severe water deficit (SWD) relative to well-watered (WW) C. canephora cv. Conilon Clone CL153 (CL153) and Coffea arabica cv. Icatu (Icatu) plants. Selected DEGs were annotated with the Gene Ontology (GO) terms: "response to water deprivation" (GO:0009414), "response to desiccation" (GO:0009269). Red: up-regulated DEGs; blue: down-regulated DEGs.

Regulation Patterns of Photosynthetic and Other Biochemical Coffee Related DEGs
Due to the fundamental role of energy pathways (photosynthesis and respiration), antioxidative protection and membrane lipid dynamics involved in coffee acclimation to environmental stresses, a specific search was also performed among the DEGs associated with these processes. Results showed that an important set of responsive DEGs (184) associated with such crucial biochemical coffee processes were affected by drought, mostly when considering SWD, where the majority of these DEGs were detected (Table S7;  . Antioxidant activity related DEGs also followed this general pattern, being mostly down-regulated in both genotypes, especially under SWD (CL153: 16, Icatu: 13). The number of these DEGs significantly increased from MWD to SWD, especially in CL153 that also increased the level of up-regulated DEGs under the harshest drought condition. By contrast, DEGs associated with cellular respiration were mostly up-regulated in CL153 under SWD (20 out of 34), while in Icatu, they were all up-regulated under MWD (all 4).

Regulation Patterns of Transcription Factors among Responsiveness DEGs
The search of Arabidopsis thaliana homologs' Transcription Factors (TFs) among the DEGs found only five TFs (Table S8). The Ethylene-responsive TF (WIN1), the TF MYB60, and the ABC transporter G family member 22 (ABCG22) were down-regulated under SWD in the two genotypes. Among the remaining TFs (all found only in CL153), the Dehydration-responsive element-binding protein 1A (DREB1A) and the NAC domaincontaining protein 55 (NAC055) were up-regulated regardless of the drought severity. In contrast, the ABC transporter G family member 22 (ABCG22) was down-regulated only under SWD.
Given that this analysis revealed a very low number of the TFs, we also searched DEGs annotated in the reference genome as TFs. With this search, a total of 83 TFs were found among DEGs, mostly in CL153 and predominantly under SWD: 26 and 62 TFs in CL153 plants, and 17 and 32 TFs in Icatu plants under MWD and SWD, respectively (Table S8). The majority of TFs were up-regulated in CL153, especially under SWD (34 up and 28 down) than under MWD (21 up and 5 down). The same pattern was reported for Icatu, with a high number of TFs being found under SWD (17 up and 15 down) and under MWD (6 up and 11 down). In both drought treatments, the Ethylene-responsive TF ERF027 followed by the Dehydration-responsive element-binding protein 1D (DREB1D) were the most up-regulated TFs in CL153 (Table S8). In Icatu, the ethylene-responsive transcription factor 15 was up-regulated under the two water deficits levels, together with the basic leucine zipper 6 under MWD and the basic leucine zipper 5 under SWD (Table S8).

Regulation Patterns of Phosphatases and Protein Kinases DEGs
Due to its importance in stress acclimation, a specific search was performed to understand the patterns of phosphatases and protein kinases among DEGs. A total of 49 phosphatase activity-DEGs were detected with a prevalence under SWD and in higher numbers in CL153 (Table S9). Approximately one-third of them were up-regulated in CL153 (13) and one-fifth in Icatu (10) under SWD, with a notable up-regulation of the phosphatase 2C 74 in CL153 and the Major allergen (Mal D1) in Icatu (Table S9).
Drought had a significant impact on 272 protein kinases (Table S9). Under MWD, a similar number of protein kinases was detected in the two genotypes (38), while SWD had a higher impact on protein kinases that were mostly down-regulated by this stress (181 in CL153 vs. 96 in Icatu). Besides an uncharacterized kinase protein, SWD triggered the most up-regulation of the G-type lectin S-receptor-like serine/threonine-protein kinase (RLK1) in CL153, which was also the most up-regulated kinase in Icatu (Table S9).

Enriched GO Terms of Drought-Related DEGs
Overall, in both genotypes, there was an increase in enriched GO terms as drought severity increased, being predominant in down-regulated DEGs ( Figure 5). Results revealed a very specific response to drought, with the catalytic activity being the only category commonly enriched in both drought treatments and both genotypes. Under MWD, GO terms categories were predominantly enriched in Icatu (16 in Icatu vs. 5 in CL153). Even so, only CL153 showed an enrichment in up-regulated DEGs linked to catalytic and DNAbinding TF activities ( Figure 5). SWD had a stronger impact on a high number of categories (30 in CL153 and 31 in Icatu), with 20 of them being commonly altered in both genotypes. Under this stress, only CL153 showed enriched categories in up-regulated DEGs that were linked to the nucleus, DNA-binding TFs, and oxidoreductase activities.

Impacts of Water Deficit on the Transcriptomic Profile of Coffea canephora (CL153) and C. arabica (Icatu)
The present study demonstrated that water deficit alters the transcriptome profile in the two coffee genotypes. Thousands of expressed genes were identified and annotated ( Figure 1; Table S1) in agreement with previous drought studies in coffee [29,52,53]. However, contrary to those studies, we examined the effects of a long-term drought experiment in which stress has been imposed gradually. In addition, this study had the advantage of comparing two genotypes from different coffee species that have been found to display distinct physiological resilience to SWD conditions, particularly at the photosynthetic functioning and apparatus components [24,25]. Our findings showed that only a small number of genes was affected by MWD in the two genotypes (CL153: 5%; Icatu: 4% DEGs), whereas a substantial impact was observed under SWD (CL153: 17%; Icatu: 10% DEGs; Table S2; Figures 2 and 3). In comparison, in the drought-sensitive C. canephora cv. Conilon 109 only 0.59% of genes responded to drought, while 3.63% were found in the drought-tolerant C. canephora cv. Conilon 120, in a 14-day drought experiment [29]. Therefore, the high responsiveness values found here under MWD (in a drought level similar to [29]) suggests a high drought tolerance in the coffee cultivars of this study (CL153 and Icatu). In fact, under MWD, the potential functioning of the photosynthetic apparatus of these plants was not significantly impaired [19], although some protective mechanisms (e.g., zeaxanthin and HSP70 content, antioxidative activity) already begun to be reinforced (Semedo et al., unpublished data).
The genes responsive to drought also differed significantly between the two coffee genotypes (Tables 1-4). In CL153, reticuline oxidase genes were predominant in the two water deficit levels, including the most up-regulation of PAP20. Several antioxidant genes such as peroxidase 4, thioredoxin, and FAD-related genes as reticuline oxidases have been previously identified in this genotype as being involved in stress acclimation, including high temperatures [34]. Thus, the up-regulation of PAP20, an acid phosphatase involved in hydrolase activity and metal ion binding, helps to alleviate the PAP stress signal that usually accumulates during drought and light stress, inducing the expression of stress-responsive genes [54], which could regulate the impact of water deficit in CL153. Reticuline oxidase genes were also up-regulated in CL153 even under MWD, together with cytochrome P450, a carotenoid cleavage dioxygenase, and the TF ERF027 (Table 1). The up-regulation of these genes, together with the TF ERF027, is probably linked to a protective mechanism of CL153 to avoid oxidative damage, as previously documented in this genotype [55]. Carotenoids are essential components of the photosynthetic apparatus, being susceptible to oxidation processes that break the carotenoid backbone [56]. This cleavage reaction is catalyzed by carotenoid cleavage dioxygenases leading to apocarotenoids that usually arise through the attack of ROS [57]. Apocarotenoids have an essential role in abiotic stress response, acting as precursors of ABA synthesis that coordinates plant responses to stress, namely stomata closure to minimize water loss during drought and therefore suggests the activation of ABA-related mechanisms to minimize drought in CL153, as shown in other drought-tolerant coffee genotypes [29]. This is in line with the increase of ABA synthesis found under MWD and SWD in CL153 (although also in Icatu) plants [19]. Specifically, the strong down-regulation of the ROP-interactive CRIB motif-containing protein 4 under SWD also suggests the involvement of ABA in the response of drought in this genotype. These types of genes are usually involved in the interaction between auxin-and ABA-regulated processes, which often show an antagonistic effect, that is, positively regulating auxin signaling while negatively regulating ABA signaling [58].
In contrast, several heat shock genes were up-regulated in Icatu even under MWD, as the glycoside hydrolase family 79, in line with the higher abundance of the heat shock protein 70 kDa (HSP70) (unpublished data). Under SWD, genes involved in binding, transporter, and transferase activities were up-regulated, namely the UDP-glycotransferase 75D1. These genes are involved in a high number of developmental processes and stress responses, including cell wall modification, plant hormone activation, and the production of antioxidants in response to stresses, including drought [59,60], which would counterbalance the down-regulation of TFs (ORG2) under the two water deficits (see also below). In coffee, as in other plants, ROS can be scavenged by several enzymes and non-enzyme antioxidants, such as ascorbic acid, glutathione, carotenoids, phenolic compounds, ascorbate peroxidase, or catalase [61], which are usually strongly reinforced in Icatu under drought [62] and would explain why this genotype can endure even the effects of extreme water deficit [25], beyond the transcriptomic results found here. Under drought, particularly under SWD, we found a down-regulation of DEGs associated with important physiologically and biochemically related processes, such as photosynthesis, lipid metabolism, cellular respiration, and antioxidant activity (except the cellular respiration in Icatu under MWD; Figure 4). The impact of the SWD level on the transcripts of the two genotypes was also revealed by the enrichment analysis that showed a prominent down-regulation of GO terms involved in general metabolic processes, integral components of the membrane, and especially in catalytic activities ( Figure 5). In fact, up-regulated enriched GO terms were only found in CL153 and involved catalytic and DNA-binding TF activities under MWD, while under SWD, the nucleus, oxidoreductase, and DNA-binding TF activities were enriched ( Figure 5). However, previous physiological, biochemical, and proteomic studies revealed an almost insensitivity of Icatu to the severe drought impact, at least in the C-assimilatory pathway functioning, as well as a stronger triggering of protective molecules [24,25,62]. This suggests the involvement of post-transcriptional regulations and indicates the need for complementary integrated studies considering several layers of plant response, from physiology and biochemistry to molecular levels. Such studies would be crucial to understand transcriptomic findings in the context of plant acclimation to environmental constraints [34].

Role of Aspartic Proteases and Protectant Proteins in Water Deprivation and Desiccation in Coffee
Among the 22 DEGs involved in response to drought, ten were related to the Aspartic Protease in Guard Cell 1 (ASPG1; Table 5). These genes have been identified in different plant species being the ASPG1 usually involved in plant acclimation to drought stress [63,64]. Over the last decade, an increasing number of publications have highlighted the involvement of aspartic proteases in plant defense responses against a diversity of abiotic and biotic stresses [65,66]. For instance, plants overexpressing aspartic proteases as APA1 have been shown to be more tolerant to water deficit due to changes in stomatal behavior induced by the regulation of the ABA signaling pathway [67]. In Arabidopsis, the ASPG1 has been shown to be involved in drought stress resistance, in addition to its role in the degradation of seed storage proteins [68,69]. Arabidopsis mutants overexpressing ASPG1 were shown to recover more efficiently from drought since ASPG1 led to a significant increase in ABA sensitivity by guard cells and activation of antioxidant enzymes that prevent plants from oxidative damage [68]. A gene homologous to ASPG1 from potato has also been shown to be down-regulated under drought and up-regulated upon re-watering, suggesting a role under drought stress [65]. Thus, the up-regulation of ASPG1 in CL153 would help this genotype to mitigate the effects of drought. By contrast, a down-regulation of this gene was observed in Icatu, even though a previous study found an increase in ABA, and a strong stomatal conductance reduction in this genotype [25], which are related to the gene expression found here.
Drought tolerance in the coffee genotypes also involved the up-regulation of Lea14 in CL153 (under the two water deficits) and in Icatu (under SWD). Late embryogenesis abundant proteins (Lea) have been found in a large number of plants being up-regulated during osmotic stress, where they bind to enzymes to prevent loss of activity functioning as cellular stabilizers during stress conditions [70]. Lea are also expressed under water deprivation conditions and associated with improved drought tolerance by inducing rapid stomatal closure [71,72]. Thus, the overexpression of Lea14 could also be involved in the water scarcity response of these two coffee genotypes.

Drought-Responsiveness Transcription Factors
In addition to protective proteins, drought also triggered a high number of TFs in the two coffee genotypes (e.g., 88 TFs,) that were mainly up-regulated (Table S8). This upregulation of TFs in response to drought was more frequent under SWD than under MWD and higher in CL153 than in Icatu, suggesting an important role in regulating drought stress tolerance, potentially through improved cellular protection, and particularly in the first genotype. TFs usually regulate genes involved in cellular protection from stress damage (e.g., osmoprotectants, antioxidants), as well as signal transduction and transcriptional regulation [73]. In coffee, several TFs have been previously observed to respond to drought exposure. For instance, a wide-throughout transcriptomic study found only 22 probable responsive TFs, namely from the Myb superfamily [29], being the overexpression of TFs, including the DREB gene family also reported [74][75][76]. Likewise, in our study, four TFs were found to be significantly up-regulated in CL153 under drought (ERF027, DREB1D, DREB1A, and NAC055 in both water deficits) and three others in Icatu (TF15 under the two water deficits, plus basic leucine zipper 6 under MWD, and zipper 5 under SWD) suggesting their involvement in the drought-response tolerance of these genotypes. In fact, ERF TFs usually regulate responses to abiotic stresses, including cold, drought, heat, and salt, being also involved in hormone signaling and hormone-mediated stress-responses through stress phytohormones as ABA [77,78] and ethylene [79]. TFs as DREBs are induced upon drought imposition, positively regulating drought-responsive genes [63,80] such as the Lea previously reported in this study. DREB/CBF belongs to the ERF (ethyleneresponsive element binding factors) superfamily of TFs that play a pivotal role in adaptation to biotic and abiotic stresses [81]. NAC TFs are also relevant in ABA and ethylene pathways responding to drought stress [82]. When cells are under water deficit, ABA accumulates and binds to soluble receptors. These results in the release and activation of Open Stomata genes (members of protein kinases) that phosphorylate basic leucine zipper TFs to control gene expression in the nucleus [83]. Therefore, the overexpression of these genes could sustain an enhanced drought tolerance [84], as observed in these two coffee genotypes [25]. Expression of ERFs can be induced by ethylene and ABA under biotic and abiotic stresses that also interact with other plant hormone pathways, such as those regulated by salicylic acid or gibberellins, and suggesting coordinated interaction of hormone signaling pathways to regulate the expression of TF genes during stress responses [81]. ERFs also seem to regulate ROS-responsive genes, resulting in decreased accumulation of ROS and enhancing tolerance to multiple abiotic stresses such as drought, salt, and freezing [85].

DEGs Involved in Phosphatases and Protein Kinases Affected by Drought
Protein phosphatases and kinases are major post-translational regulators of numerous cellular processes and signaling networks [86]. Protein kinases pathways are activated by sequential phosphorylation leading to the regulation of TFs and protective enzymes in response to several stresses, including drought [87]. A high up-regulation of phosphatase genes was reported under SWD (Table S9), particularly in CL153, suggesting the involvement of these genes in the drought response. The phosphatase 2C 74 that was overexpressed in CL153 under SWD is part of a major group of protein phosphatases in plants, having important roles in various biological processes [88]. Several studies have shown that PP2 genes are involved in the regulation of the ABA signaling pathway by modulating kinase activities in response to abiotic stresses [89]. Notably, the up-regulation of the major allergen Mal D1 phosphatase in Icatu under SWD (Table S8) could raise a significant concern since this gene was initially thought to be a major allergen in several fruits [90]. However, other studies showed that these proteins are also synthesized in response to biotic and abiotic stresses [91,92] and, thus, it would be interesting to determine if resistance to stresses could have consequences in terms of the allergenicity of the agronomic product.
G-type lectin S-receptor-like serine/threonine protein kinases are positive regulators of plant tolerance to several stresses [93,94]. For instance, transcriptomic analyses in foxtail millet have shown their participation in the drought tolerance response [95]. Thus, the highest up-regulation of RLK1 in both CL153 and Icatu under SWD agrees with its role in drought regulation in coffee, most likely also triggered by ABA signaling pathways induced by protein phosphatases as referred above. Ultimately, the primary response beyond ABA levels will be the stomatal closure in coffee leaves, regulating ion transport in guard cells and decreasing drought severity in coffee genotypes as reported in other species [96][97][98].

Coping with Drought: Lessons from Crossed Transcriptomic, Physiological, and Biochemical Studies
Climate changes are expected to include an increased frequency of water scarcity events, including the intensity of severity and duration, posing a growing threat to coffee's global supply chain. Agroforestry systems can be useful to mitigate some of these harsh effects, reducing the risk of coffee production losses and contributing to the sustainability of crops [18,31]. However, to maintain the global supply of coffee, it is also important to promote the screening and development of tolerant varieties to face the increasingly expected impacts of droughts.
The novel mRNA-Seq study reported here provided information on qualitative and quantitative differences between the two cultivars, CL153 and Icatu. Overall, comparative transcriptome analysis led to the identification of drought-responsive genes and genotypedependent genes responsible for the different drought tolerance responses, the TFs ERF027, DREB1D, and the basic leucine zipper genes, as well as genes linked to water deprivation and desiccation, the ASPG1 and the protectant Lea14. These genes will be essential for future crop improvement programs, such as the development of drought-resilient coffee varieties. However, some of the transcriptomic results found here do not fully agree with earlier physiological and biochemical studies showing a greater tolerance of Icatu than CL153 under SWD. Icatu has previously been shown to be barely affected by drought, showing minor impacts on photosynthetic functioning (e.g., Amax, Fv/Fm) and components (e.g., electron carriers, photosystems, and ribulose-1,5-bisphosphate carboxylase oxygenase (RuBisCO) activity) under SWD, in contrast to CL153 which was clearly affected. This better performance of Icatu under such harsh water shortage conditions was associated with the triggering of mechanisms of (thermal) energy dissipation and ROS control over the photosynthetic machinery [24,25,62]. Such enrichment of detox pathways was also accompanied by metabolic and proteomic changes in Icatu, which included the reinforcement of thylakoid electron transport rates and some electron carriers, and the triggering of protective cyclic electron transport involving both photosystems. These processes would help maintain the photochemical use of energy while controlling the presence of reactive molecules of chlorophyll and oxygen [24,25,107]. In this context, the existence of strong post-transcriptional regulations is very likely, probably involving alternative splicing, noncoding RNAs (including siRNA, miRNA, lncRNA), RNA-binding proteins (RBPs) as well as protein modifications [108]. In the future, combining these technological methods with bioinformatic tools and physiological experiments will allow a more holistic insight into the regulation and control of biological processes in coffee.

•
Even though drought had an impact on the leaf transcriptome of both coffee genotypes, our results revealed that both genotypes are more drought-resistant than other coffee cultivars. • Drought triggered a specific response associated with the magnitude of water deficit, which was also genotype-dependent since few DEGs and pathways were common to treatments and both genotypes. By comparison, MWD only had a minor effect on the transcripts of both genotypes.

•
There was a predominance of protective genes (more in CL153) associated with antioxidant activities, including genes involved in water deprivation and desiccation, such as Lea and aspartic proteases. • A significant number of TFs, including ERF, DREB, and the leucine zipper, were found to be significantly up-regulated under drought. Together with the large number of phosphatases and protein kinases we found, these results suggest the involvement of ABA signaling in the drought tolerance of these genotypes.

•
Coupled with the previous physiological and metabolic results, our study provides novel and timely information showing several layers of response and suggesting the existence of post-transcriptional regulations in the two coffee genotypes, which should be further investigated.