Transcriptomic Leaf Profiling Reveals Differential Responses of the Two Most Traded Coffee Species to Elevated [CO2]

Abstract

As atmospheric [CO₂] continues to rise to unprecedented levels, understanding its impact on plants is imperative to improve crop performance and sustainability under future climate conditions. In this context, transcriptional changes promoted by elevated CO₂ (eCO₂) were studied in genotypes from the two major traded coffee species: the allopolyploid Coffea arabica (Icatu) and its diploid parent, C. canephora (CL153). While Icatu expressed more genes than CL153, a higher number of differentially expressed genes were found in CL153 as a response to eCO₂. Although many genes were found to be commonly expressed by the two genotypes under eCO₂, unique genes and pathways differed between them, with CL153 showing more enriched GO terms and metabolic pathways than Icatu. Divergent functional categories and significantly enriched pathways were found in these genotypes, which altogether supports contrasting responses to eCO₂. A considerable number of genes linked to coffee physiological and biochemical responses were found to be affected by eCO₂ with the significant upregulation of photosynthetic, antioxidant, and lipidic genes. This supports the absence of photosynthesis down-regulation and, therefore, the maintenance of increased photosynthetic potential promoted by eCO₂ in these coffee genotypes.

1. Introduction

Coffee is among the most important crops worldwide, being grown in the tropical region in approximately 80 countries where it plays a determinant economic and social role [1]. The livelihoods of about 25 million farmers’, mostly smallholders, depend on this highly labor-intensive crop that involves 100–125 million people in its worldwide chain of value [2,3]. Trading is dominated by two species: the allotetraploid Coffea arabica L. (2n = 4x = 44) and one of its diploid ancestors, C. canephora Pierre ex A. Froehner, which together are responsible for nearly 99% of the global coffee production [1,2,3]. Coffea arabica appears to have originated from a single polyploidization that occurred in very recent evolutionary times (< 50,000 years ago) in the plateaus of Central Ethiopia, between two diploid species: C. eugenioides and C. canephora [4,5]. The two parental species are closely related, and the two subgenomes have a very low sequence divergence (i.e., 1.3% average difference for genes [6]). Coffea arabica displays disomic inheritance with bivalent pairing of homoeologous chromosomes [4,5], but its ability to tolerate a broader temperature range than its diploid parents is not related to the overall higher expression of the homoeologous genes [7]. Moreover, the extremely low levels of genetic variation found in C. arabica [8] when compared to its diploid progenitors is a concern in the context of environmental changes and measures are needed to increase the environmental, economic, and social sustainability of coffee.

The continuous rapid rise in ambient air [CO₂] (aCO₂) raises many questions related to food security, namely how this will impact the quality and yield of crops [9,10]. Elevated air [CO₂] (eCO₂) is implicated in climate changes through warming and altered precipitation patterns, which in turn are associated to higher evapotranspiration, soil salinization, pest and disease pressure, and reduced availability of arable soil [11,12,13,14]. Still, due to a C-fertilization effect that might increase crop yield and plant tolerance to environmental stresses, positive eCO₂ impacts can occur through an increase of photosynthesis [1,15,16,17,18,19].

Despite the potential implications of eCO₂, the impacts at physiological, biochemical, and structural levels on the two main producing coffee species were only recently studied [13,20,21]. Leaf physiological studies revealed some differences between genotypes, but no clear overarching pattern of species-dependent responses to eCO₂ occurs. Overall, the exposure to eCO₂ supports important net photosynthesis increases between 34% and 49% under environmental-controlled conditions at 700 μmol mol⁻¹ [20], ca. 40–60% at 550 μmol mol⁻¹ under field conditions using a Free Air CO₂ Enrichment (FACE) facility [13], or even above 60% on Open Top Chambers (OTC) systems [22], under natural or semi-controlled environmental conditions, respectively. In fact, under eCO₂, no photosynthetic and stomatal down-regulation is detected (the so-called negative acclimation) in coffee plants in contrast to a wide range of plant species, where photoassimilate accumulation and lower investment in photosynthetic components can compromise the stimulation effect of eCO₂ on photosynthesis [13,20]. Additionally, eCO₂ improves the physiological status of coffee plants, which allows a greater coffee tolerance to heat [23,24] and drought [22], contributing to preserve coffee bean quality under supra-optimal temperature conditions [24].

Identifying the genes underlining plant responses to altered environmental conditions is clearly crucial to better understand crop sustainability in a scenario of climate changes [25]. In particular, the transcriptomic approach is a powerful tool to detect genes and pathways involved in stress tolerance [26]. In other species, transcriptomic studies revealed an alteration in the patterns of gene expression when plants are exposed to eCO₂, affecting key processes, such as photosynthesis [27,28], and respiration [29]. Very recently, transcriptomic studies also identified gene expression differences when contrasting populations were exposed to ambient and eCO₂ [30,31], showing acclimation responses to this important changing environmental condition [32]. Molecular studies to understand Coffea species responses to climate changes have been fostered by the sequencing of the genome of C. canephora, yielding 25,574 annotated protein-coding genes [33]. In this context, this work aimed to elucidate the transcriptomic responses of Coffea leaves to eCO₂ (700 μmol mol⁻¹), which is forecast to occur in the near future, considering two important cropped cultivars in Brazil, from the two main producing species: C. arabica cv. Icatu (Arabica type of coffee) and C. canephora cv. Conilon Clone 153 (Robusta type of coffee). This transcriptomic approach complements the thorough mineral, physiological, and biochemical leaf analysis performed under the same CO₂ conditions studied here [18,20,23], allowing a broader and integrated picture regarding the mechanisms behind Coffea responses, and a new glimpse on the future of this crop under changing environmental conditions.

2. Results 2.1. Overall Transcriptome Profiling and Mapping Statistics

Deep RNA sequencing of Coffea generated an average of 27.8 million raw reads (Table S1). After stringent quality assessment and data filtering, we obtained 93% of clean reads corresponding to an average of 25.9 million clean reads (Table S1). Overall, a high proportion of uniquely reads was aligned to the reference genome of Coffea canephora (unmapped reads: CL153: 2.8%; Icatu: 3.2%), demonstrating a very high coverage over the species transcriptome. In Icatu, the proportion of mapped reads aligned to a unique position was similar between aCO₂ and eCO₂, whereas in CL153, this proportion increased under eCO₂ (77.7% vs. 83.9%). After checking for replicate variance in the expression profile with PC analysis, the outlier replicate 1C was removed from further analysis (Figure S1).

2.2. Differential Gene Expression in Response to eCO₂

Icatu expressed 21,714 genes under aCO₂ and 21,659 genes under eCO₂, which is higher than the number of genes expressed by CL153 in those conditions (aCO₂: 20728) or at eCO₂ (21186) (Figure 1; Table S2). Nevertheless, as a response to eCO₂, Icatu expressed 4895 annotated DEGs (23%) while a higher number was found in CL153 (6486; 31%) considering a normalized non-zero log2 fold change expression and a false discovery rate FDR < 0.01. Annotated DEGs were slightly up-regulated in Icatu (52%) and slightly down-regulated in CL153 (51%) (Figure 1). Change values of DEGs ranged from −8.08 to 6.11 in Icatu (Table S3) and from −5.31 to 9.88 to CL153 (Table S4). It is worth nothing the expression of 24% of DEGs in Icatu and 26% of DEGs in CL153, which retrieved no annotation using the Uniprot genes and the local database (Table S2).

A high number of DEGs (2799) were found to be commonly expressed by the two genotypes as a response to eCO₂ (Figure 2), corresponding to 57% of all DEGs in Icatu, and to 43% of DEGs in CL153.

2.3. Major CO₂-Responsive DEGs in the Two Genotypes

In Icatu, eCO₂ led to the overexpression of 12 DEGs with fold-changes higher than four (FC > 4; Table 1). They were involved in the acid metabolic process (Cc01_g09180), unsaturated fatty acid biosynthetic process (Cc01_g05170), photosystem II repair, proteolysis, and thylakoid membrane organization (Cc07_g12510), metal binding (Cc05_g12830), cell redox homeostasis (Cc08_g16300), chloroplast organization (Cc07_g04750), circadian rhythm, regulation of transcription, and response to stress (Cc07_g12820), and two aquaporins involved in water transport (Cc02_g05510 and Cc04_g17080), plus three uncharacterized proteins (Table 1).

In the opposite regulation, 10 DEGs showed an FC < -4 under eCO₂, being involved in DNA-binding transcription factor activity (Cc02_g03420), protein ubiquitation (Cc06_g10420, Cc01_g14750), carbohydrate catabolic process (Cc02_g27160, Cc_g07710), and electron transfer (Cc06_g08260, Cc09_g06040), plus three uncharacterized proteins (Table 1).

In CL153, 19 DEGs with FC > 4 were expressed under eCO₂ (Table 2). Those up-regulated DEGs were mainly involved in general biosynthetic processes (Cc05_g05400, Cc02_g25790, Cc09_g03820, Cc05_g13060), as well as in general defense responses to stresses (Cc02_g03420, Cc03_g12230, Cc06_g13950, Cc03_g15270, Cc10_g12850), cell wall organization (Cc01_g08230), glyoxylate/tricarboxylic acid cycle (Cc05_g04940), and oxidoreductase activity (Cc01_g09500), plus seven uncharacterized ones. In comparison, eCO₂ led to a strong down-regulation of four DEGs with FC < 4 that were involved in phytoalexins biosynthesis (Cc11_g04610) and transmembrane transporter activity (Cc07_g08450), plus two uncharacterized proteins (Table 2).

2.4. Differential Gene Expression Responses between Icatu and CL153

Icatu expressed a higher number of genes than CL153 at aCO₂ (more 986 genes) or at eCO₂ (more 473 genes; Figure 1; Table S2). Nevertheless, there was always a common response by the two genotypes, either at eCO₂ or aCO₂: 3451 genes (1704 down- and 1747 up-regulated) were commonly expressed between genotypes under both aCO2 and eCO2 (Figure 3). These corresponded to 51% of all DEGs found at eCO₂ and to 59% of all DEGs found at aCO₂, while the remaining DEGs were exclusive to only one CO₂ level. Fold change values between genotypes varied from -13.36 to 14.23 under aCO₂ (Table S5) and from -10.99 to 12.17 under eCO₂ (Table S6).

However, the two genotypes responded differently at the same air [CO₂] level. Under aCO₂, CL153 expressed 20 up-regulated DEGs with FC > 10 involved in general processes, such as the lipid transport (Cc00_g10190), root and shoot development (Cc00_g08140), chitin catabolic process (Cc05_g00740), flavonoid metabolic process (Cc00_g25490), cell cycle (Cc03_g11280), UDP-glucose metabolic process (Cc00_g21940), and transferase activity (Cc00_g07170; Table 3). Interestingly, more than half (12/20) are uncharacterized proteins.

In contrast, the most down-regulated DEGs in CL153, or in other words the most up-regulated ones in Icatu, were involved in general stress responses (Cc00_g30240, Cc00_g17430, Cc00_g20380, Cc08_g08870, Cc07_g14720), flavonoid metabolic process (Cc00_g29880), and oxidoreductase activity (Cc09_g05150, Cc05_g05110) plus eight uncharacterized proteins (Table 3).

Under eCO₂, when compared to Icatu, CL153, most DEGs were involved in lipid transport (Cc00_g10190), leaf morphogenesis and shoot formation (Cc08_g09540), cellular differentiation (Cc01_g04210), and response to oxidative stress (Cc11_g08120), plus three uncharacterized ones (Table 4). In contrast, DEGs in Icatu were involved in oxidoreductase activity (Cc05_g05110), and alkaloid biosynthetic process (Cc00_g11770) but mainly in response to stresses (Cc00_g17430, Cc00_g30240, Cc00_g06270), plus three uncharacterized ones (Table 4).

2.5. Common DEGs to Both Genotypes but with Opposite Patterns of Expression to eCO₂

Of the DEGs commonly expressed by the two genotypes under eCO₂, 83% (2799) showed a similar pattern of regulation while the remaining 17% (567) showed opposite patterns (Figure 4). From those, 311 DEGs were down-regulated in CL153 but up-regulated in Icatu while 255 showed the opposite pattern (Figure 4). The highest change was found in Cc00_g06550, a gene involved in the gibberellin biosynthetic process and leaf development that was down-regulated in CL153 but up-regulated in Icatu and in Cc02_g03430, a dehydration-responsive element that was up-regulated in CL153 but down-regulated in Icatu (Table S7).

2.6. Significantly Enriched GO Categories of the Two Genotypes in Response to eCO₂

Functional characterization of eCO₂-responsive DEGs in Icatu found 3923, 4097, and 4168 gene hits in biological process, molecular function, and cellular components GO terms, respectively (Table S8). However, only three GO terms, all from the biological process category, were significantly enriched in Icatu under eCO₂ (FDR < 0.05, p < 0.001), having the highest normalized enrichment score (Figure 5; Table S9). Each enriched GO term was associated with 18 to 36 up-regulated DEGs. In comparison, the functional characterization of CL153 found 5205, 5373, and 5590 of CO₂-responsive gene hits in biological process, molecular function, and cellular components GO terms, respectively, which is higher than the ones reported for Icatu.

Enriched GO terms were more significantly enriched in CL153 than in Icatu, with 10 of them (6 from biological process and 4 from molecular function GO terms) being significantly up-regulated and 3 (2 from biological process and 1 from cellular component GO terms) significantly down-regulated in CL153 under eCO₂ (Figure 5; Table S10). Each enriched GO term was associated with 5 to 65 up-regulated DEGs and with 11 to 46 down-regulated DEGs, respectively.

Two categories of GO terms were significantly affected by eCO₂ but with opposite effects on the two genotypes. The Molecular Function GO:0016679 (oxidoreductase activity acting on diphenols and related substances as donors) was significantly enriched among up-regulated DEGs in CL153 relative to Icatu under eCO₂. The Cellular Component GO:0042170 (plastid membrane) was significantly enriched among down-regulated DEGs in CL153 relative to Icatu under eCO₂ (Figure 6).

2.7. Significant Metabolic Pathways

In Icatu, several metabolic pathways were found to be enriched in up-regulated and down-regulated genes as a response to eCO₂, although none was significant, either considering the KEGG or WikiPathways database (Table 5). In contrast, two KEGG pathways were found to be significantly enriched in CL153, both up-regulated: carotenoid biosynthesis and cutin, suberine, and wax biosynthesis (Table 5).

WikiPathways found six enriched pathways in CL153 although none was significant. No metabolic pathways were found to be significantly enriched either at aCO₂ or eCO₂ when comparing the two genotypes with the KEGG database (Figure S2). However, WikiPathways found a significantly enriched pathway at aCO₂—genetic interactions between sugar and hormone signaling -, detected as being up-regulated in CL153 in comparison to Icatu under aCO₂ (Table 5).

2.8. Identification of Photosynthetic and Other Biochemical-Related Responsive DEGs

As a response to eCO₂, a total of 171 and 190 DEGs were found associated respectively with Icatu and CL153 responses to environmental modifications at photosynthetic and other related biochemical components, although responses differ between the two genotypes (Table S11). Photosynthesis and chlorophyll metabolic process-related DEGs were mostly down-regulated under eCO₂ in CL153, while Icatu only showed a slight down-regulation of these DEGs (CL153: 78% and 70%; Icatu: 51% and 53%; Figure 7). GO-terms associated with RuBisCO were largely up-regulated in both genotypes (both 67%) under eCO₂ (Table S11).

Antioxidant components were slightly up-regulated under eCO₂ and with similar levels found in both genotypes (CL153: 57%; Icatu: 51%). Lipid metabolic FAD and LOX genes also showed a similar trend of upregulation in both genotypes (CL153: 67%; Icatu: 63%), as well as the GO terms associated to cellular respiration and pyruvate kinase that were up-regulated, especially in Icatu (CL153: 69% and 75%; Icatu: 82% and 100%, respectively). However, malate dehydrogenase-related DEGs were more up-regulated in CL153 than in Icatu (CL153: 83%; Icatu: 50%). 3. Discussion

Transcriptome responses of plants to eCO₂ vary according to the species and even to the genotypes studied. For instance, in wheat (Triticum aestivum cv. Norstar), 1022 genes were differentially expressed in response to eCO₂ (700 μmol mol⁻¹; [34]). In contrast, at 700 μmol mol⁻¹, Plantago lanceolata expressed only 33-131 DEGs, although numbers increased to 689-853 in populations adapted to eCO₂ [32]. These transcriptionally activated genes are thought to help plants to cope with eCO₂, although the underlying genomic mechanisms that are affected, modulated, and controlled by eCO₂ are still not fully understood [35].

In this study, we found that Icatu expressed more genes than CL153 either at aCO₂ or eCO₂ (Figure 1), with values similar to the ones previously reported for C. canephora [36] and C. arabica [37] under aCO₂. Beyond well-established roles in promoting new physiological, metabolic, and morphological traits [38], polyploidy can promote non-uniform genome and transcriptome alterations [7,39], which might explain the differences found here between the polyploid C. arabica and the diploid, C. canephora and supports a differential response by these species to environmental conditions.

3.1. Transcriptomic Responses of Icatu to eCO₂

Significant changes were detected in DEGs linked to increased plant tolerance and adaptation to abiotic factors in Icatu (Table 1). For instance, the glycotransferase, Cc01_g15810, up-regulated in this genotype under eCO₂, belongs to a large family of genes involved in plant growth and response to biotic and abiotic stresses [40,41]. Glycotransferases catalyze the transfer of sugar moieties from active sugar molecules to a variety of acceptor molecules, namely lipids, which allow plant cells to modulate their biochemical properties [40]. The Cc01_g0517 also up-regulated in Icatu is an integral membrane protein that enhances the content of dienoic fatty acids, which in turn might increase tolerance to abiotic factors [42]. In accordance with the upregulation of these genes, changes in the lipid profile of Icatu were reported to increase under eCO₂, which is quite relevant since qualitative and quantitative lipid adjustments, namely of chloroplast membranes, were identified as crucial to the ability shown by Icatu for long-term acclimation to heat and high responsiveness to eCO₂ [43]. Other DEGs that were up-regulated under eCO₂ included transcriptional factors involved in chloroplast organization (Cc07_g04750) and in the repair of PSII, as well as in the maintenance of the thylakoid structure (Cc07_g12510). The expression of these genes is frequently related to the prevention of photoinhibition and irreversible damage of PSII that limit the growth and development of plants [44,45]. Indeed, a high preservation was found on the performance of the photosynthetic apparatus of Icatu at high temperatures, under eCO₂ rather than under aCO₂ [18]. Aquaporins, members of major intrinsic proteins (MIPs), were also found to be up-regulated in Icatu (Cc02_g05510 and Cc04_g17080; Table 1). Among other roles, aquaporins control the hydraulic conductance of cells and organs, managing water flux in and out of veins and stomatal guard cells, being part of the response to drought in plants, namely in coffee leaves [22]. They can also increase membrane permeability to CO₂ in mesophyll and stomatal guard cells, the latter increasing the effectiveness of RuBisCO, potentially influencing transpiration efficiency [46]. In coffee leaves, aquaporin transcriptional genes seem not to significantly change the regulation of stomatal conductance (g_s), since it is barely affected under well-watered conditions and eCO₂ [13,20]. However, in another C. arabica genotype, faster stomatal closure rates (in well-watered conditions) and the maintenance of plant hydraulic conductance (in drought conditions) observed under eCO₂ were linked to a higher transcript abundance of most aquaporin genes [22], which can be seen as a preventive advantage if coffee plants are afterwards exposed to drought conditions.

It is also worth mentioning the down-regulation of protein ubiquitination genes (Cc01_g14750, Cc06_g10420) in Icatu, as well as genes involved in electron transfer (Cc09_g06040, Cc06_g08260) and DNA binding (Cc02_g03420), which are usually involved in stress signaling, launching, and repairing damage components [47]. In fact, an increase of these genes is usually expected for removing misfolded or damaged proteins that accumulate as an exposure to abiotic stress [48]. Therefore, its down-regulation is consistent with the stress relief impact of eCO₂ found in Coffea species regarding drought [22] and heat [18,23] stresses, likely associated with the greater photochemical use of energy due to a higher carboxylation rate [13,20,22]. This would also decrease the photorespiration rate and the oxidative pressure at the chloroplast level, reducing the reactive oxygen species (ROS) formation and the excitation level in the photosystems [49,50].

In accordance, several genomic, proteomic, and physiological studies under eCO₂ showed an increase in metabolic processes linked to carbohydrate availability, cellular growth, photosynthetic carbon fixation, and stress defense in C3 plants [51], which would explain the three GO terms significantly enriched in Icatu under eCO₂ (secondary metabolic processes, plastid organization, and response to light intensity). Even though no specific pathways were found to be significantly enriched in Icatu as a response to eCO₂ (Table 5), a global increase was observed regarding these processes. Notably, these alterations on the gene expression profile of Icatu under eCO₂ are not reflected on bean physical and chemical characteristics, which was previously found to show marginally, if at all, changes in response to eCO₂ [24]. This was also suggested for caffeine in other C. arabica genotypes, using leaf caffeine contents as a proxy for concentrations in the beans [52].

3.2. Transcriptomic Responses of CL153 to eCO₂

In response to eCO₂, DEGs in CL153 were mainly involved in general biosynthetic processes (Cc05_g05400, Cc02_g25790, Cc09_g03820, Cc05_g13060) and responses to stresses (Cc06_g13950, Cc03_g15270, Cc10_g12850; Table 2) that have been shown to protect plants against general abiotic stresses as heat and ozone damage [53,54,55]. For example, the overexpression of Isoprene Synthase (ISPS), Cc05_g05400 and Cc05_g04940 homologs gene, enhances the production of leaf chlorophyll and carotenoid contents and has a significant effect on leaf and inflorescence growth in Arabidopsis [56] although eCO₂ might inhibit isoprene emissions in other plants [57]. The overexpression of the aldehyde decarbonylase CER3 (Cc02_g25790) under eCO₂ could also be an important adaptation due to its role in the synthesis of major wax components [58], influencing stomatal development [59] and providing extra energy for photosynthesis [60], which, as discussed before, showed no down-regulation in coffee leaves under eCO₂. Another example is the ethylene-responsive transcription factor ERF025 (Cc02_g03420, Cc03_g12230), an important biotic and abiotic stress-responsive gene [61,62] involved in transcriptional regulatory mechanisms to eCO₂ [56]. The upregulation of cytochrome P450s (Cc01_g09500) and the glucomannan 4-beta-mannosyltransferase (Cc01_g08230) gene also reflects important transcriptomic changes in CL153 under eCO₂ as they play important roles in plant acclimation/tolerance against several abiotic stresses, promoting cell wall organization and increasing its expansion [63]. In contrast, the Cc07_g08450 gene involved in transmembrane transporter activity and the Momilactone A synthase gene (Cc11_g04610) were both down-regulated under eCO₂ (Table 2). The opposite pattern is usually found as a response to abiotic and biotic stresses since these genes are thought to be involved in plant tolerance [64,65,66].

According to the overexpression of these genes in CL153, significant enriched GO terms were related to general biological processes as cell wall biogenesis, secondary metabolic processes, cell wall macromolecule metabolic processes, plant-type cell wall organization or biogenesis, and polysaccharide metabolic processes, as well as molecular functions linked with oxidoreductase and transferase activity (Figure 5). Previous studies showed that eCO₂ alters plant structure by inducing changes in the rate of cell division and expansion, and in the accumulation of starch within chloroplasts to accomplish a higher photosynthetic rate, which are helped by a high activity of oxireductase and transferase genes [67,68]. Indeed, some components of the photosynthetic machinery in CL153 are up-regulated under eCO₂, with starch usually increasing 69% under eCO₂ and showing the highest thylakoid electron transport capability among tested coffee genotypes [20].

Not surprisingly, the cutin, suberine, and wax KEGG pathway was found to be significantly enriched in CL153 under eCO₂ (Table 5) due to the overexpression of the genes mentioned above. Indeed, plants under eCO₂ may generate more biomass and height than those under aCO₂ [3,69], playing an important ecological role as a form of carbon sequestration. The carotenoid KEGG biosynthesis pathway was also significantly enriched in CL153 (Table 5), which is also explained by the extra carbohydrates promoted by eCO₂ that provide both the elemental substrate and energy source to support such synthesis in C3 plants as the coffee tree [70]. In fact, at ambient temperature (25/20 °C), eCO₂ promotes higher contents of several carotenoids in CL153, being particularly significantly for neoxanthin and carotenes [23].

3.3. Main Transcriptomic Differences between the Two Genotypes

Even though the two genotypes are closely related due to the allopolyploid origin of C. arabica involving C. canephora as a progenitor [71,72], significant transcriptomic differences were found between the two genotypes in response to eCO₂, and even at aCO₂ (Figure 3), with relevant fold change differences found between them (aCO₂: -13.36 to 14.23; eCO₂: -10.99 to 12.17). In fact, despite the close relationship, C. canephora and C. arabica evolution and selection determined different ecological requirements (namely regarding, altitude and temperature, as well as rainfall amount), and distinct acclimation responses to environmental stresses (e.g., to drought, heat, cold) (for a review, see [73,74]), which would be associated to different transcriptional profiles, as it is the case in the present work. Icatu always expressed more genes than CL153 at aCO₂ or eCO₂, which were involved in different metabolic processes than those of CL153. At aCO₂, there was an overexpression of DEGs involved in general metabolic processes (e.g., lipid transport, leaf and growth development, chitin catabolic process, flavonoid and UDP-glucose metabolic process), while Icatu showed an overexpression of DEGs associated with general stress responses (Table 3). Under eCO₂, CL153 DEGs were mostly involved in general growth and cell differentiation processes, while Icatu DEGs were mainly stress related (Table 4).

Almost half of DEGs were commonly expressed by the two genotypes under eCO₂, with the majority (83%) showing a similar pattern of regulation while 17% showed opposite patterns of expression (Table S7). The higher upregulation of Cc00_g06550 in Icatu in relation to CL153 could indicate an important role to boost development and growth of Icatu under eCO₂ since gibberellins are a class of hormones well known to regulate plant growth and development and increasing tolerance to abiotic stresses [75]. In contrast, the upregulation of the dehydration-responsive element Cc02_g03430 in CL153 and its down-regulation in Icatu as a response to eCO₂ could indicate a potentially better tolerance response of CL153 than Icatu, since this gene family includes key transcription factors involved in plant responses to various abiotic stresses, including drought tolerance in coffee leaves [76,77].

No metabolic KEGG pathways were found to be significantly enriched when comparing the two genotypes, regardless of [CO₂] (Figure S2). However, WikiPathways found a significantly enriched pathway (genetic interactions between sugar and hormone signaling) among up-regulated DEGs in CL153 in comparison to Icatu at aCO₂, which is involved in the hub-regulation of plant growth, development, and survival [78]. The absence of significant differences is of great importance to interpret biological results regarding the acclimation of the coffee plant to eCO₂, suggesting that this new environmental condition maintained the balance of most functional pathways, and could help to overcome environmental stresses, such as high heat and drought.

3.4. Transcriptomic Responses Associated to Physiological and Biochemical Behavior of Coffee

Several studies have shown that coffee trees respond positively to eCO₂ through increased leaf photosynthetic rates, carbohydrate, lipids, and even increased yield under adequate water supply [18,23,43]. For this reason, particular attention was given to GO terms related to photosynthesis, respiratory, lipid metabolism, and antioxidative pathways that resulted in the identification of 171 DEGs in Icatu and 190 DEGs in CL153. Photosynthesis and chlorophyll metabolic process-related DEGs were mostly down-regulated in CL153, but only marginally down-regulated in Icatu under eCO₂ (Figure 7). Despite this, GO terms associated with RuBisCO were largely up-regulated in Icatu and in CL153 under eCO₂. This supports physiological results showing increased total carboxylase RuBisCO activity (confirmed both through enzyme activity and A-to-internal CO₂ curves) under eCO₂ in Icatu, together with the moderate reinforcement of other photosynthetic components (e.g., associated to thylakoid electron transport) and in the photosynthetic capacity, A_max [20]. This reflects an absence of photosynthetic negative acclimation in coffee plants grown under eCO₂ [18,20,22], in sharp contrast to many other species [16,35]. Contrasting results from biochemical and gene expression analysis were also found in these genotypes submitted to high temperatures (42ºC), in which maximal gene expression above 8- (Icatu) and 250-fold (CL153) were concomitant with minimal ascorbate peroxidase activity [23]. Therefore, our findings further highlight the need of complementary studies using physiological to molecular tools to evaluate coffee responses to eCO₂.

The promotion to greater photo-assimilate availability under elevated CO₂ usually leads to an upregulation of the genes associated with the respiratory pathway [29,79], and the presence of respiratory glycolysis intermediates may enhance the respiration potential [51]. This is in accordance with the upregulation found in cellular respiration-, pyruvate kinase-, and malate dehydrogenase-related DEGs (Figure 7).

This study also found an upregulation of DEGs linked to antioxidant components (Figure 7). In accordance, some antioxidative enzymes (Cu, Zn-superoxide dismutase, and ascorbate peroxidase) and other protective molecules (heat shock protein 70, HSP70) have been shown to be reinforced in these genotypes under eCO₂ [23]. This defensive and protective role of antioxidant components has been reported to reduce the oxidative stress under eCO₂ also in other plants [80,81,82]. Notably, a reduction in antioxidative molecules under eCO2 was reported in several plant species, which is seen as reflecting a weaker ability to cope with sudden stress events by the plants grown under eCO₂ [49,83]. Instead, in coffee, these components were somewhat reinforced, together with a greater use of energy through photochemical processes, and the concomitant reduction of photorespiration by eCO₂. Altogether, these would reduce the excitation pressure in the photosystems and the reactive oxygen species (ROS) formation [49,50], as observed, namely, in Icatu [23].

Additionally, as referred above, it is recognized that coffee plant acclimation to environmental drifts, namely heat [18], high irradiance [84], and eCO₂ [43], is intimately related to the ability to perform qualitative and quantitative adjustments in the lipid matrix of chloroplast membranes. These include an altered amount in fatty acids and unsaturation degree, as well as changes in the balance between the main lipid classes. It has been shown that both genotypes, but especially Icatu, have a high lipid responsiveness under eCO₂, with a substantial increase in galactolipid classes, and low reductions in phospholipids, particularly phosphatidylglycerol, which play important roles in chloroplast membranes and, therefore, in photosynthetic functioning [43]. Notably, both genotypes presented an upregulation of FAD and LOX genes closely associated with lipid metabolism (Figure 7).

3.5. Adapting Coffee Crop for eCO₂: Links to Future Research

Photosynthesis in the coffee crop is stimulated by eCO₂, as in almost all C3 species, at adequate levels of water and temperature, but eCO₂ can also mitigate the negative impacts of warming [18] and drought [22]. Our study reports for the first time a distinct functional transcriptomic response between the two main coffee-producing species in response to eCO₂ conditions. Clear differences were found, being photosynthetic and chlorophyll metabolic genes are more affected in CL153 than in Icatu, even though the two genotypes seem to respond to eCO₂ with an increase in the activity of genes and pathways linked to acclimation / defense capabilities, as well as cell modification genes under heat [23] and drought [22,85] conditions. Therefore, it is of paramount importance to further test these and other coffee genotypes, in multiple environmental conditions expected to occur under field conditions in the near future. This would be important to strengthen the in situ management, towards a climate change mitigation that can grant a sustainable development of the coffee crop, which is already being affected by climate changes [25]. It is also noteworthy the identification of relevant expression patterns of several ‘novel’ genes coding uncharacterized/hypothetical proteins. The assignment of function to those uncharacterized genes is still experimentally a substantial challenge [86,87], but together with better coffee genomic databases related to gene expression, protein and metabolite presence, and the integration with physiological and biochemical information already available, it would be critically important to deciphering coffee responses to future environmental conditions.

3.6. Conclusions eCO₂ caused significant changes in the transcriptomic profile of the two genotypes, with DEGs being much more abundant in CL153 than in Icatu. Approximately half of DEGs were found to be a common response by the two genotypes to eCO₂ with the majority of these showing the same type of regulation between genotypes. Significant fold change differences were found between the two genotypes, being much greater than those found between [CO₂] conditions within each genotype.

The functional characterization unveiled important differences in the expression of genes between the two genotypes even though there was an absence of large significant responses showing that eCO₂ maintained the balance of most functions in coffee. Under eCO₂, CL153 mostly expressed DEGs associated to general biological processes and to a lower extent to abiotic stress responses. In contrast, Icatu increased the upregulation of genes linked to plant tolerance and adaptation to abiotic factors, such as oxidative stress control genes or aquaporins. A group of genes closely involved with adjustments in the membrane lipid matrix, in chloroplast/thylakoid organization, and PSII repair were also up-regulated in Icatu. This likely supports a greater photosynthetic performance under stressful conditions of heat [18] and drought [22], in the plants grown under eCO₂ than under aCO₂.

Altogether, our findings support previous complementary studies at physiological and biochemical levels, providing a complete picture regarding coffee (and other crops) response and sustainability under eCO2 and associated climate changes. 4. Materials and Methods 4.1. Plant Material and Experimental Design

Seedlings from the two main producing species C. canephora Pierre ex A. Froehner cv. Conilon Clone 153 (CL153) and C. arabica L. cv. Icatu Vermelho (Icatu), the latter an introgressed variety resulting from a cross between C. canephora and C. arabica cv. Bourbon Vermelho that was further crossed to C. arabica cv. Mundo Novo, were grown in 12-L pots until 1.5 years of age, in a greenhouse, under ambient [CO₂] (aCO₂). Afterwards, plants were transferred into walk-in growth chambers (EHHF 10000, ARALAB, Portugal), and grown for another 10 months in 28-L pots under controlled environmental conditions of temperature (25/20 °C, day/night), relative humidity (70–75%), irradiance (ca. 700 μmol m⁻² s⁻¹), photoperiod (12 h), and air [CO₂] of 380 μmol mol^-1 (aCO₂) or 700 μmol mol^-1 (eCO₂) [18,20]. The eCO₂ was chosen considering values that are estimated to be reached along the second half of the current century [9,10]. The plants were grown in an optimized substrate consisting of a mixture of soil, peat, and sand (3:1:3, v/v/v), and fertilized as previously described [20].

Newly matured leaves from plagiotropic and orthotropic branches were collected under photosynthetic steady-state conditions (after ca. 2h of illumination) from the upper (illuminated) part of each plant (6 per treatment), flash frozen in liquid nitrogen, and stored at −80 °C until analysis. Along the experiment, the plants were maintained without restrictions of water, nutrients, and root development, the latter judged by visual examination at the end of the experiment [20].

4.2. Total RNA, Poly(A) RNA Isolation, and Library Preparation Total RNA from the two genotypes (CL153, Icatu) × the two [CO₂] levels (aCO₂ and eCO₂) × 3 biological replicates (each one a pool from 3 different plants) was isolated using the RNeasy Plant Mini Kit (Qiagen, Germany) according to the manufacturer’s protocol and digested with DNase I using an on-column Qiagen DNase set to remove putative genomic DNA, as per the manufacturer’s instructions. The concentration and quality of each RNA sample was determined in 1.5% agarose–TBE gel electrophoresis containing GelRed Nucleic Acid Gel Stain (Biotium, CA, USA) by evaluating the integrity of the 28S and 18S ribosomal RNA bands and absence of smears. All RNA samples were individually analyzed for the possible presence of DNA contamination by standard PCR reactions (35 cycles) using primers designed for the ubiquitin (UBQ) gene, in the absence of cDNA synthesis. Total RNA concentration and purity were further verified through BioDrop Cuvette (BioDrop, UK). Samples were further quantified in an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) at Instituto Gulbenkian de Ciências, Oeiras, Portugal to verify total RNA quality. 4.3. Processing and Mapping of Illumina Reads

The messenger RNA (mRNA) libraries were constructed with the Illumina “TruSeq Stranded mRNA Sample Preparation kit” (Illumina, San Diego, CA, USA) and sequenced separately on a Hiseq 2000 (Illumina) at the MGX platform (Montpellier GenomiX, Evry Cedex, France). Raw reads were deposited in the NCBI Sequence Read Archive, BioProject accession PRJNA606444. Raw reads obtained by sequencing were quality checked using FastQC version 0.11.8 [88]. The reads were screened for contaminants using FastQ Screen version 0.13 [89] against the genome of 14 default pre-indexed FastQ putative contaminant species and adapters for both PhiX genome and TruSeq Illumina adapters. After trimming with Trimmomatic version 0.38 [90], qualified reads were mapped to the reference genome of C. canephora downloaded from the Coffee Genome Hub (http://coffee-genome.org, accessed on 13 November 2019; [33]) using STAR version 2.6.1 [91] with default settings as quantMode set to GeneCounts, as well as for producing the mapping statistic. HTSeq-count version 0.11.0 [92] was used to count how many reads mapped uniquely to each gene, discarding reads in multiple alignments, and avoiding the increase of false positives. Relevant parameters used included the default mode “union” and option “stranded = reverse”. Samtools version 1.9 [93] and gffread version 0.9.9 [94] were used throughout the analysis to convert files and obtain general statistics of the genome mapping. For exploratory analysis, Principal Coordinate Analysis (PCoA) was conducted to verify the relatedness of every replicate using the function prcomp in R software version 3.5.1 [95]. Afterwards, log2 transformation of Fragments Per Kilobase Of Exon Per Million Fragments Mapped (FPKM+1) and quantile normalization were performed using Cufflinks version 2.2.1 [96].

4.4. Identification of Differentially Expressed Genes

High-throughput RNA-sequencing combinations were performed to detect the effect of eCO₂ (aCO₂ vs. eCO₂; aCO₂ as the control) in CL153 and in Icatu. We also studied transcriptomic changes in the two genotypes at aCO₂ and at eCO₂, detecting changes in CL153 relative to Icatu (as the control). The DESeq2 version 3.8 [97] was employed to identify differentially expressed genes (DEGs). DESeq2 fits generalized linear models to each gene, using shrinkage estimation for dispersions and fold changes to improve the stability and interpretability of estimates [87]. The resulting values were adjusted using the Benjamini and Hochberg’s approach for controlling the FDR [98]. Genes with a normalized non-zero log2 fold change expression and an FDR < 0.01 were defined as differentially expressed. Venn diagrams were developed to identify common genes and specific DEGs to the effect of eCO₂ and the response of each genotype using VennDiagram v1.6.20 R package [99].

4.5. Functional Annotation and Enrichment Analysis

Functional annotations for the protein-coding genes in the C. canephora genome were downloaded from the Coffee Genome Hub (http://coffee-genome.org, accessed on 13 November 2019). Because a very high number of reads were mapped to the C. canephora genome, it was further used as the reference genome for the analyses. BLAST2GO version 1.4.4 [100] was used for mapping and functional annotation of DEGs (parameters: E-Value-Hit-Filter 1.0E−6, Annotation Cutoff 55, gene ontology [GO] Weight 5, Hsp-Hit Coverage Cutoff 20). After mapping, DEGs were filtered to exclude multiple hits to the same gene, keeping only the one that showed the highest identity percentage. The genes were characterized using GO terms of molecular function (MF), biological process (BP), and cellular component (CC). A local BLAST database was built to map DEGs to the highest identity Uniprot gene hits, using the C. canephora genome. GO mapping and GO annotation was then performed with Blast2GO command line interface using the Uniprot genes and the local database. Then, to obtain a more complete list of terms, all GO ancestors of each identified term were retrieved using QuickGO API. To visualize the size and the overlap of the gene sets described by a given GO term, a gene set enrichment analysis (GSEA) was performed using the default parameters in WebGestalt (WEB-based Gene SeT AnaLysis Toolkit) with ranked lists of DEGs, sorted by log2FC, against the Arabidopsis thaliana functional database (http://www.webgestalt.org/, accessed on 30th March 2020). A GSEA was chosen to take into account the strength of feature differentiation among DEGs, in an attempt to detect relevant changes to genes in the dataset that cause alterations in given terms or pathways from the knowledge base. Gene sets with an FDR < 0.05 were considered enriched, where the FDR is the estimated probability that a gene set with a given normalized enrichment score represents a false positive finding. Pathway enrichment analysis was also performed on WebGestalt using the KEGG and Wikipathway databases. WebGestalt results were plotted using the R ggplot2 [101] version 3.3.2 library.

To relate the transcriptomic answer of the two coffee genotypes with their physiological and biochemical responses, DEGs annotated with the GO terms referenced in [20] (e.g., photosynthesis, chlorophyll metabolic process, ribulose-1,5-bisphosphate carboxylase/oxygenase, antioxidant activity, cellular respiration, malate dehydrogenase activity, and pyruvate kinase activity) were searched using the QuickGO API [102], identifying all descendant GO terms annotated with Blast2GO. We also specifically searched for all FAD- and LOX-related proteins, which have been reported as the most important for lipid profile dynamics related to stress acclimation in coffee plants [43].

Figure 1. Effect of CO2 enhancement in Icatu and CL153. Left: Number of total expressed genes at aCO2 (yellow) and at eCO2 (blue) and number of annotated differentially expressed genes (DEGs; green) at eCO2 vs. aCO2 considering a normalized non-zero log2 fold change expression and an FDR < 0.01, in Icatu (above) and CL153 (below). Right: Number of annotated up- and down-regulated significant DEGs at eCO2 vs. aCO2 in Icatu (above) and CL153 (below).

Figure 2. Down- (left) and up-regulated (right) annotated differentially expressed genes (DEGs) commonly expressed by the two genotypes at eCO2 (green), specific of CL153 (yellow) and specific of Icatu (blue).

Figure 3. Down- (left) and up-regulated (right) annotated differentially expressed genes (DEGs) in CL153 vs. Icatu. Numbers indicate specific genotype differences in DEGs at eCO2 (yellow), at aCO2 (pink) and shared between both [CO2] levels (brown).

Figure 4. Expression patterns of DEGs commonly expressed by both genotypes to eCO2. Pie chart: = indicates DEGs with the same pattern of regulation between genotypes while ≠ indicates DEGs with opposite regulations between genotypes. Bars indicate DEGs up-regulated in one genotype and down-regulated in the other.

Figure 5. Gene Set Enrichment Analysis (GSEA) performed by WEB-based Gene SeT AnaLysis Toolkit (WebGestalt). Significantly (FDR < 0.05, p< 0.001) enriched Gene Ontology (GO) terms among down- and up-regulated differentially expressed genes (DEGs) ranked by increasing log2 fold-change (FC), considering the effect of eCO2 in CL153 (left) and Icatu (right). Gene Ontology (GO) terms are grouped by the main category - Biological Process (GO:BP), Molecular Function (GO:MF), and Cellular Component (GO:CC). Counts indicate the number of DEGs annotated with each term/pathway and dots are colored by ascending normalized enrichment score (NES).

Figure 6. Gene Set Enrichment Analysis (GSEA) performed by WEB-based Gene SeT AnaLysis Toolkit (WebGestalt). Significantly (FDR < 0.05, p< 0.001) enriched Gene Ontology (GO) terms among down- and up-regulated differentially expressed genes (DEGs) ranked by increasing log2 fold-change (FC), considering the response between CL153 and Icatu genotypes under aCO2 (left) and eCO2 (right). Gene Ontology (GO) terms are grouped by main category - Molecular Function (GO:MF) and Cellular Component (GO:CC). Counts indicate the number of DEGs annotated with each term/pathway and dots are colored by ascending normalized enrichment score (NES).

Figure 7. Proportion of significantly up- (blue) and down-regulated (yellow) DEGs associated to physiological and biochemical responses in coffee that were found at eCO2 vs. aCO2, in Icatu and CL153. DEGs were searched in the following GO terms: Photosynthesis, Chlorophyll Metabolic Process, RuBisCO, Antioxidant Activity, FAD and LOX (Lipid Metabolic Process), Cellular Respiration, Malate Dehydrogenase (MDH) activity, and Pyruvate Kinase (PK) activity.

Gene	Protein Name	FC	Biological Processes
Top up-regulated
Cc01_g09180	Flavonol 7-O-beta-glucosyltransferase	6.11	acid metabolic process
Cc10_g08360	Uncharacterized protein	5.60
Cc01_g05170	Omega-6 fatty acid endoplasmic reticulum isozyme 2	5.60	unsaturated fatty acid biosynthetic process
Cc07_g12510	ATP-dependent zinc metalloprotease	5.19	photosystem II repair, proteolysis, thylakoid membrane organization
Cc02_g07390	Uncharacterized protein	5.18
Cc05_g12830	Zinc finger Constans-Like 10	4.86	metal binding
Cc08_g16300	Uncharacterized protein At5g39865	4.83	cell redox homeostasis
Cc07_g04750	Zinc finger Constans-Like 2	4.63	chloroplast organization
Cc07_g12820	LHY	4.42	circadian rhythm, regulation of transcription, response to stress
Cc02_g05510	Aquaporin PIP subfamily	4.36	water transport
Cc04_g17080	Aquaporin, MIP family, PIP subfamily	4.07	water transport
Cc09_g02890	Uncharacterized protein	4.05
Top down-regulated
Cc02_g03420	Ethylene-responsive transcription factor ERF025	−5.11	DNA-binding transcription factor activity
Cc02_g03180	Uncharacterized protein	−5.21
Cc06_g10420	E3 ubiquitin- ligase CIP8	−5.46	protein ubiquitination
Cc11_g07700	Uncharacterized protein	−5.62
Cc02_g27160	Vicianin hydrolase (Fragment)	−5.69	carbohydrate catabolic process
Cc11_g07710	Two-component response regulator	−5.77	carbohydrate catabolic process
Cc06_g08260	Umecyanin	−5.83	electron transfer activity, metal ion binding
Cc09_g06040	Pentatricopeptide repeat-containing At2g35130	−6.30	electron transfer activity, metal ion binding
Cc11_g12580	Uncharacterized protein	−7.69
Cc01_g14750	Chloroplast import apparatus 2	−8.08	protein ubiquitation

Gene	Protein Name	FC	Biological Process
Top up-regulated
Cc05_g05400	Isoprene chloroplastic	9.88	isoprene synthase activity
Cc05_g04940	Isocitrate lyase	8.60	glyoxylate cycle, tricarboxylic acid cycle
Cc01_g05170	Uncharacterized protein	8.46
Cc05_g00220	Uncharacterized protein	8.26
Cc02_g25790	Very-long-chain aldehyde decarbonylase CER3	8.00	alkane biosynthetic process
Cc02_g20300	Uncharacterized protein	7.83
Cc02_g03420	Ethylene-responsive transcription factor ERF025	7.79	defense response to fungus, ethylene-activated signaling pathway
Cc03_g12230	Desiccation-related PCC13-62	7.59	cellular response to desiccation, cellular response to salt
Cc06_g13950	EG45-like domain containing 1	7.56	cellular response to hypoxia
Cc03_g15270	Desiccation-related PCC13-62	7.55	cellular response to desiccation, cellular response to salt
Cc01_g08230	Glucomannan 4-beta-mannosyltransferase 9	6.96	cell wall organization
Cc01_g09500	Cytochrome P450 94A1	6.89	oxidoreductase activity
Cc09_g03820	12-oxophytodienoate reductase 1	6.88	oxylipin biosynthetic process
Cc02_g08150	Uncharacterized protein	6.54
Cc00_g12860	Uncharacterized protein	6.40
Cc10_g12850	Stearoyl	6.13	cellular response to hypoxia
Cc00_g11660	Uncharacterized protein	6.12
Cc05_g13060	(-)-germacrene D synthase	6.11	terpenoid biosynthetic process
Cc07_g14440	Uncharacterized protein	6.08
Top down-regulated
Cc11_g04610	Momilactone A synthase	-5.01	phytoalexins biosynthesis
Cc06_g10420	Uncharacterized protein	-5.09
Cc11_g07700	Uncharacterized protein	-5.22
Cc07_g08450	Detoxification 6	-5.31	transmembrane transporter activity

Gene	Protein Name	FC	Biological Processes
Top up-regulated
Cc11_g10870	Uncharacterized protein	14.23
Cc00_g10190	Acyl- -binding domain-containing 4	12.51	lipid transport, response to ethylene, response to jasmonic acid, response to light stimulus
Cc10_g14870	Uncharacterized protein	12.03
Cc00_g08140	Homeobox-leucine zipper protein	11.89	root, shoot development, response to light, red, far-red light phototransduction
Cc05_g00740	Acidic endochitinase	11.87	chitin catabolic process, polysaccharide catabolic process
Cc00_g25490	8-hydroxyquercetin 8-O-methyltransferase	11.57	flavonoid metabolic process, methylation
Cc03_g11280	Cell division control 48 homolog B	11.29	cell cycle, cell division, protein transport
Cc07_g01500	Uncharacterized protein	11.05
Cc01_g01970	F-box LRR-repeat At4g14096	11.00
Cc00_g17790	Uncharacterized protein	10.92
Cc00_g26610	Uncharacterized protein	10.87
Cc10_g14880	Uncharacterized protein	10.82
Cc09_g08250	Uncharacterized protein	10.71
Cc03_g09260	Uncharacterized protein	10.59
Cc00_g21940	UTP--glucose-1-phosphate uridylyltransferase	10.54	UDP-glucose metabolic process
Cc01_g01360	Uncharacterized protein	10.54
Cc06_g16070	Uncharacterized protein	10.36
Cc11_g00660	Uncharacterized protein	10.24
Cc09_g08930	Uncharacterized protein	10.22
Cc00_g07170	UDP-glucose flavonoid 3-O-glucosyltransferase 3	10.00	transferase activity
Top down-regulated
Cc07_g14720	1-aminocyclopropane-1-carboxylate oxidase 2	−10.06	stress response, ethylene biosynthetic process
Cc09_g05150	Cytochrome P450 89A2	−10.09	oxidoreductase activity
Cc00_g29880	8-hydroxyquercetin 8-O-methyltransferase	−10.14	flavonoid metabolic process, methylation
Cc08_g08870	BURP domain RD22	−10.25	response to salt stress
Cc11_g04040	Uncharacterized protein	−10.28
Cc10_g04410	Uncharacterized protein	−10.40
Cc00_g23750	Uncharacterized protein	−10.48
Cc10_g14350	Uncharacterized protein	−10.66
Cc00_g20380	FK506-binding 2	−10.70	stress response
Cc05_g05110	Cytochrome P450 71D10	−10.70	oxidoreductase activity
Cc00_g17430	Isoflavone reductase homolog P3	−10.76	response to cadmium ion, response to oxidative stress
Cc06_g14770	Uncharacterized protein	−10.97
Cc11_g00400	Uncharacterized protein	−11.08
Cc00_g29670	Uncharacterized protein	−11.54
Cc04_g13840	Uncharacterized protein	−11.60
Cc00_g30240	Cysteinease inhibitor 5	−13.36	cellular response to heat, defense response

Gene	Protein Name	FC	Biological Processes
Top up-regulated
Cc00_g10190	Acyl- -binding domain-containing 4	12.17	lipid transport, response to ethylene, response to light stimulus
Cc10_g14870	Uncharacterized protein	11.60
Cc04_g16250	Uncharacterized protein	11.30
Cc08_g09540	Carotenoid cleavage dioxygenase chloroplastic	11.21	leaf morphogenesis, secondary shoot formation
Cc05_g00740	Uncharacterized protein	10.67
Cc11_g08120	Plant cadmium resistance 2	10.62	response to oxidative stress
Cc01_g04210	Homeobox-leucine zipper ATHB-8	10.57	cell differentiation
Top down-regulated
Cc00_g06270	Annexin D8	−10.07	response to cold, response to heat, response to salt stress, response to water deprivation
Cc00_g30240	Cysteinease inhibitor 5	−10.09	cellular response to heat, defense response
Cc00_g17430	Isoflavone reductase homolog P3	−10.23	response to oxidative stress
Cc02_g38970	Uncharacterized protein	−10.28
Cc00_g19610	Uncharacterized protein	−10.39
Cc00_g11770	Tabersonine 16-O-methyltransferase	−10.44	alkaloid biosynthetic process
Cc00_g01740	Uncharacterized protein	−10.77
Cc05_g05110	Cytochrome P450 71D10	−10.99	oxidoreductase activity

Database	ID	Description	Counts	NES	p-value	FDR
eCO2 vs. aCO2 | CL153
KEGG	map00906	carotenoid biosynthesis	12	1.99	<0.001	8.12 × 10−4
	map00073	cutin, suberine and wax biosynthesis	5	1.76	1.69 × 10−3	4.91 × 10−2
CL153 vs. Icatu | eCO2
WikiPathways	WP3661	genetic interactions between sugar and hormone signaling	12	1.67	2.15 × 10−2	2.90 × 10−2

Supplementary Materials

Supplementary materials can be found at https://www.mdpi.com/1422-0067/21/23/9211/s1. Figure S1. Principal component (PC) analysis of the raw reads generated by RNA-sequencing with (A) and without replicate 1C (B). Figure S2. Enriched KEGG metabolic pathways found among DEGs in CL153 vs. Icatu at aCO2 (A) and eCO2 (B) showing the highest (blue) and lowest (yellow) normalized enrichment scores. Dark colors indicate statistically significant KEGG metabolic pathways under FDR <0.05. Table S1. Genome mapping showing the alignment and reads counting results of the transcriptome of Icatu and CL153 against the genome of Coffee canephora. A, B, C correspond to the biological replicates. RAW READS: number of reads obtained after sequencing. CLEAN READS: number of reads passing the Illumina quality filters and downstream filters. % CLEAN: percentage of reads passing filters compared to the number of raw reads. UNIQUE: number of reads aligned to a unique position. % UNIQUE: proportion of reads aligned to a unique position compared to the number of clean reads. MULTIPLE MAP: number of reads aligned to exons of several overlapping genes. % MULTIPLE MAP: proportion of reads aligned to exons of several overlapping genes compared to the number of clean reads. % UNMAPPED: proportion of non-aligning reads compared to the number of clean reads. Table S2. Total gene expression and differentially expressed genes (DEGs) of Icatu and CL153 samples, grown at aCO2 and aCO2. Genotype: Number of total expressed genes at aCO2 and at eCO2 in Icatu and CL153 genotypes, number of all, unknown and annotated DEGs at eCO2 vs. aCO2, with the respective number of up- and down-regulated DEGs. [CO2]: Number of total expressed genes by Icatu and by CL153 in aCO2 and eCO2, number of all, unknown and annotated and annotated DEGs in CL153 vs. Icatu, with the respective number of up- and down-regulated DEGs. Table S3. List of down- and up-regulated Icatu DEGs at eCO2 vs. aCO2 with gene identification, protein name and fold-change (FC) values. Genes are sorted by fold-change in descending order. Table S4. List of down- and up-regulated CL153 DEGs at eCO2 vs. aCO2 with gene identification, protein name and fold-change (FC) values. Genes are sorted by fold-change in descending order. Table S5. List of down- and up-regulated CL153 DEGs when comparing to Icatu, at aCO2 with gene identification, protein name and fold-change (FC) values. Genes are sorted by fold-change in descending order. Table S6. List of down- and up-regulated CL153 DEGs when comparing to Icatu, at eCO2 with gene identification, protein name and fold-change (FC) values. Genes are sorted by fold-change in descending order. Table S7. List of commonly regulated DEGs at eCO2 vs. aCO2 in both CL153 and Icatu, but with opposite patterns of expression. Table S8. Lists of all differentially expressed genes (DEGs) in Icatu and CL153 genotypes under aCO2 and eCO2 that were used as input for the Gene Set Enrichment Analysis (GSEA). Uniprot IDs represent the top blastx hit of each gene. Log2FC (FC) values were calculated for the following comparisons: eCO2 vs. aCO2 in Icatu (Icatu_700vs380), eCO2 vs. aCO2 in CL153 (CL153_700vs380), CL153 vs. Icatu under aCO2 (380_CL153vsIcatu) and CL153 vs. Icatu under eCO2 (700_CL153vsIcatu). Table S9. Gene Set Enrichment Analysis (GSEA) of differentially expressed genes (DEGs), considering the effect of eCO2 in Icatu and CL153, performed with WEB-based Gene SeT AnaLysis Toolkit (WebGestalt). Significantly enriched Gene Ontology (GO) terms from each category - Biological Process (BP), Molecular Function (MF) and Cellular Component (CC). Values indicate the number of DEGs annotated with each term and pathway (Counts), normalized enrichment scores (NES), p-value and False Discovery Rate (FDR <0.05; p <0.001). Table S10. Gene Set Enrichment Analysis (GSEA) of differentially expressed genes (DEGs), considering the response between CL153 and Icatu genotypes under aCO2 and eCO2. Significantly enriched Gene Ontology (GO) terms from each category -Molecular Function (MF) and Cellular Component (CC). Values indicate the number of DEGs annotated with each term and pathway (Counts), normalized enrichment scores (NES), p-value and False Discovery Rate (FDR <0.05; p <0.001). Table S11. Icatu and CL153 photosynthetic-related responsive DEGs to eCO2. The searched GO-terms included: "Photosynthesis", "Chlorophyll metabolic process", "Ribulose-bisphosphate carboxylase activity", "Antioxidant activity", "Lipid metabolic process (LOX, FAD)", "Cellular respiration", "Malate dehydrogenase activity", and "Pyruvate kinase activity". Gene and protein annotation, main function(s) retrieved from UniprotKB, and respective fold-change. Blue indicates up-regulation and yellow indicates down-regulation.

Author Contributions

Conceptualization, A.I.R.-B., and J.C.R.; Data curation, I.M., I.F., P.H.D. and J.C.R.; Formal analysis, I.M., I.F. and P.H.D.; Funding acquisition, J.C.R. and A.I.R.-B.; Investigation, L.F.G., A.S.F., F.C.L., J.C.R. and A.I.R.-B.; Methodology, I.R.-B. and J.C.R.; Project administration, J.C.R. and A.I.R.-B.; Supervision, I.M., O.S.P., J.C.R. and A.I.R.-B.; Writing-original draft, I.M., I.F., P.H.D., O.S.P., F.C.L., F.M.D., J.C.R. and A.I.R.-B.; Writing-review & editing, I.M., I.F., P.D., O.S.P., F.M.D., J.C.R. and A.I.R.-B. All authors have read and agreed to the published version of the manuscript.

Funding

This work received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 727934 (project BreedCAFS), and from national funds from Fundação para a Ciência e a Tecnologia (FCT), Portugal, through the project PTDC/ASP-AGR/31257/2017, and the research units UIDB/00239/2020 (CEF), UIDP/04035/2020 (GeoBioTec), UID/AGR/04129/2020 (LEAF), and UIDB/04551/2020 (Green-it). Fellowships from the Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brazil (CNPq), and the Fundação de Amparo à Pesquisa do Estado de Minas Gerais, Brazil (FAPEMIG, project CRA-RED-00053-16), to F.M. DaMatta, are also greatly acknowledged.

Acknowledgments

The authors would like to thank Novadelta - Comércio e Indústria de Cafés Lda., as well as Paula Alves for technical assistance.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. DaMatta, F.M.; Avila, R.T.; Cardoso, A.A.; Martins, S.C.V.; Ramalho, J.C. Physiological and agronomic performance of the coffee crop in the context of climate change and global warming: A review. J. Agric. Food Chem. 2018, 66, 5264-5274.

2. Osorio, N. The Global Coffee Crisis: A Threat to Sustainable Development; International Coffee Organization: London, UK, 2002; p. 4.

3. Semedo, J.N.; Rodrigues, W.P.; Dubberstein, D.; Martins, M.Q.; Martins, L.D.; Pais, I.P.; Rodrigues, A.P.; Leitão, A.E.; Partelli, F.L.; Campostrini, E.; et al. Coffee Responses to Drought, Warming and High [CO2] in a Context of Future Climate Change Scenarios. In Theory and Practice of Climate Adaptation; Climate Change Management Series; Alves, F., Leal, W., Azeiteiro, U., Eds.; Springer: Berlin/Heidelberg, Germany, 2018; pp. 465-477. ISBN 978-3-319-72874-2.

4. Lashermes, P.; Combes, M.C.; Robert, J. Molecular characterisation and origin of the Coffea arabica L. Genome. Mol. Gen. Genet. 1999, 261, 259-266.

5. Tesfaye, K.; Borsch, T.; Govers, K.; Bekele, E. Characterization of Coffea chloroplast microsatellites and evidence for the recent divergence of C. arabica and C. eugenioides chloroplast genomes. Genome 2007, 50, 1112-1129.

6. Cenci, A.; Combes, M.C.; Lashermes, P. Genome evolution in diploid and tetraploid Coffea species as revealed by comparative analysis of orthologous genome segments. Plant Mol. Biol. 2012, 78, 135-145.

7. Combes, M.-C.; Dereeper, A.; Severac, D.; Bertrand, B.; Lashermes, P. Contribution of subgenomes to the transcriptome and their intertwined regulation in the allopolyploid Coffea arabica grown at contrasted temperatures. New Phytol. 2013, 200, 251-260.

8. Scalabrin, S.; Toniutti, L.; Di Gaspero, G.; Scaglione, D.; Magris, G.; Vidotto, M.; Pinosio, S.; Cattonaro, F.; Magni, F.; Jurman, I.; et al. A single polyploidization event at the origin of the tetraploid genome of Coffea arabica is responsible for the extremely low genetic variation in wild and cultivated germplasm. Sci. Rep. 2020, 10, 4642.

9. Stocker, T.F.; Qin, D.; Plattner, G.-K.; Tignor, M.M.B.; Allen, S.K.; Boschung, J.; Nauels, A.; Xia, Y.; Bex, V.; Midgley, P.M. Climate Change 2013. Physical Science Basis; Summary for Policymakers, Technical Summary and Frequent Asked Questions. Part of the Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; IPCC: Geneva, Switzerland, 2013; p. 203. ISBN 978-92-9169-138-8.

10. IPCC. Climate Change 2014. In Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Edenhofer, O., Pichs-Madruga, R., Sokona, Y., Farahani, E., Kadner, S., Seyboth, K., Adler, A., Baum, I., Brunner, S., Eickemeier, P., et al., Eds.; Cambridge University Press: Cambridge, UK; New York, NY, USA, 2014; p. 1435. ISBN 978-1-107-65481-5.

11. Tubiello, F.N.; Soussana, J.-F.; Howden, M. Crop and pasture response to climate change. Proc. Natl. Acad. Sci. USA 2007, 104, 19686-19690.

12. Wassmann, R.; Nelson, G.C.; Peng, S.B.; Sumfleth, K.; Jagadish, S.V.K.; Hosen, Y.; Rosegrant, M.W. Rice and global climate change. In Rice in the Global Economy: Strategic Research and Policy Issues for Food Security; Pandey, S., Byerlee, D., Dawe, D., Dobermann, A., Mohanty, S., Rozelle, S., Hardy, B., Eds.; International Rice Research Institute: Manila, Philippines, 2010; pp. 411-432.

13. Ghini, R.; Torre-Neto, A.; Dentzien, A.F.M.; Guerreiro-Filho, O.; Iost, R.; Patrício, F.R.A.; Prado, J.S.M.; Thomaziello, R.A.; Bettiol, W.; DaMatta, F.M. Coffee growth, pest and yield responses to free air CO2 enrichment. Clim. Chang. 2015, 132, 307-320.

14. Magrach, A.; Ghazoul, J. Climate and pest-driven geographic shifts in global coffee production: Implications for forest cover, biodiversity and carbon storage. PLoS ONE 2015, 10, e0133071.

15. Polley, H.W. Implications of atmospheric and climate change for crop yield. Crop Sci. 2002, 42, 131-140.

16. Long, S.P.; Ainsworth, E.A.; Rogers, A.; Ort, D.R. Rising atmospheric carbon dioxide: Plants FACE the future. Annu. Rev. Plant Biol. 2004, 55, 591-628.

17. Kirschbaum, M.U.F. Does enhanced photosynthesis enhance growth? Lessons learned from CO2 enrichment studies. Plant Physiol. 2011, 155, 117-124.

18. Rodrigues, W.P.; Martins, M.Q.; Fortunato, A.S.; Rodrigues, A.P.; Semedo, J.N.; Simoes-Costa, M.C.; Pais, I.P.; Leitao, A.E.; Colwell, F.; Goulao, L.; et al. Long-term elevated air [CO2] strengthens photosynthetic functioning and mitigates the impact of supra-optimal temperatures in tropical Coffea arabica and C. canephora species. Glob. Chang. Biol. 2016, 22, 415-431.

19. Wang, A.; Lam, S.K.; Hao, X.; Li, F.Y.; Zong, Y.; Wang, H.; Li, P. Elevated CO2 reduces the adverse effects of drought stress on a high-yielding soybean (Glycine max (L.) Merr.) cultivar by increasing water use efficiency. Plant Physiol. Biochem. 2018, 132, 660-665.

20. Ramalho, J.C.; Rodrigues, A.P.; Semedo, J.N.; Pais, I.P.; Martins, L.D.; Simões-Costa, M.C.; Leitã, A.E.; Fortunato, A.S.; Batista-Santos, P.; Palos, I.M.; et al. Sustained photosynthetic performance of Coffea spp. under long-term enhanced [CO2]. PLoS ONE 2013, 8, e82712.

21. Rakocevic, M.; Ribeiro, R.V.; Marchiori, P.E.R.; Filizola, H.F.; Batista, E.R. Structural and functional changes in coffee trees after 4 years under free air CO2 enrichment. Ann. Bot. 2018, 121, 1065-1078.

22. Avila, R.T.; Cardoso, A.; de Almeida, L.; Costa, L.C.; Machado, K.L.G.; Barbosa, M.L.; de Souza, R.P.B.; Oliveira, L.A.; Batista, D.S.; Martins, S.C.V.; et al. Coffee plants respond to drought and elevated [CO2] through changes in stomatal function, plant hydraulic conductance, and aquaporin expression. Environ. Exp. Bot. 2020, 177, 104148.

23. Martins, M.Q.; Rodrigues, W.P.; Fortunato, A.S.; Leitao, A.E.; Rodrigues, A.P.; Pais, I.P.; Martins, L.D.; Silva, M.J.; Reboredo, F.H.; Partelli, F.L.; et al. Protective response mechanisms to heat stress in interaction with high [CO2] conditions in Coffea spp. Front. Plant Sci. 2016, 7, 947.

24. Ramalho, J.C.; Pais, I.P.; Leitão, A.E.; Guerra, M.; Reboredo, F.H.; Máguas, C.; Carvalho, M.L.; Scotti-Campos, P.; Ribeiro-Barros, A.I.; Lidon, F.C.; et al. Can elevated air [CO2] conditions mitigate the predicted warming impact on the quality of coffee bean? Front. Plant Sci. 2018, 9, 287.

25. Bunn, C.; Läderach, P.; Rivera, O.O.; Kirschke, D.A. Bitter cup: Climate change profile of global production of Arabica and Robusta coffee. Clim. Chang. 2015, 129, 89-101.

26. Tallis, M.J.; Lin, Y.; Rogers, A.; Zhang, J.; Street, N.R.; Miglietta, F.; Karnosky, D.F.; Angelis, P.D.; Calfapietra, C.; Taylor, G. The transcriptome of Populus in elevated CO2 reveals increased anthocyanin biosynthesis during delayed autumnal senescence. New Phytol. 2010, 186, 415-428.

27. Taylor, G.; Street, N.R.; Tricker, P.J.; Sjodin, A.; Graham, L.; Skogstrom, O.; Calfapietra, C.; Scarascia-Mugnozza, G.; Jansson, S. The transcriptome of Populus in elevated CO2. New Phytol. 2005, 167, 143-154.

28. De Souza, A.P.; Gasper, M.; Da Silva, E.A. Elevated CO2 increases photosynthesis, biomass and productivity, and modifies gene expression in sugarcane. Plant Cell Environ. 2008, 31, 1116-1127.

29. Leakey, A.D.B.; Ainsworth, E.A.; Bernacchi, C.J.; Rogers, A.; Long, S.P.; Ort, D.R. Elevated CO2 effects on plant carbon, nitrogen, and water relations: Six important lessons from FACE. J. Exp. Bot. 2009, 60, 2859-2876.

30. Grossman, J.D.; Rice, K.J. Contemporary evolution of an invasion grass in response to elevated atmospheric CO2 at a Mojave desert FACE site. Ecol. Lett. 2014, 17, 710-716.

31. Horgan-Kobelski, T.; Matesanz, S.; Sultan, S.E. Limits to future adaptation in the invasive plant Polygonum cespitosum: Expression of functional and fitness traits at elevated CO2. J. Hered. 2016, 107, 1-9.

32. Watson-Lazowski, A.; Lin, Y.; Miglietta, F.; Edwards, R.J.; Chapman, M.A.; Taylor, G. Plant adaptation or acclimation to rising CO2? Insight from first multigenerational RNA?Seq transcriptome. Glob. Chang. Biol. 2016, 22, 3760-3773.

33. Denoeud, F.; Carretero-Paulet, L.; Dereeper, A.; Droc, G.; Guyot, R.; Pietrella, M.; Zheng, C.; Alberti, A.; Anthony, F.; Aprea, G.; et al. The coffee genome provides insight into the convergent evolution of caffeine biosynthesis. Science 2014, 345, 1181.

34. Kane, K.; Dahal, K.P.; Badawi, M.A.; Houde, M.; Hüner, N.P.A.; Sarhan, F. Long-term growth under elevated CO2 suppresses biotic stress genes in non-acclimated, but not cold-acclimated winter wheat. Plant Cell Physiol. 2013, 54, 1751-1768.

35. Xu, Z.; Jiang, Y.; Zhou, G. Response and adaptation of photosynthesis, respiration, and antioxidant systems to elevated CO2 with environmental stress in plants. Front. Plant Sci. 2015, 6, 701.

36. Guedes, F.A.; Nobres, P.; Rodrigues, D.C.F.; Menezes-Silva, P.E.; Ribeiro-Alves, M.; Correa, R.L.; DaMatta, F.M.; Alves-Ferreira, M. Transcriptional memory contributes to drought tolerance in coffee (Coffea canephora) plants. Environ. Exp. Bot. 2018, 147, 220-233.

37. Ivamoto, S.T.; Reis, O.; Júnior Domingues, D.S.; dos Santos, T.B.; de Oliveira, F.F.; Pot, D.; Leroy, T.; Esteves Vieira, L.G.; Carazzolle, M.F.; Pereira, G.A.G.; et al. Transcriptome analysis of leaves, flowers and fruits perisperm of Coffea arabica L. reveals the differential expression of genes involved in raffinose Biosynthesis. PLoS ONE 2017, 12, e0169595.

38. Iannicelli, J.; Guariniello, J.; Tossi, V.E.; Regalado, J.J.; Di Ciaccio, L.; van Baren, C.M.; Álvarez, S.I.P.; Escandón, A.S. The "polyploid effect" in the breeding of aromatic and medicinal species. Sci. Hortic. 2020, 260, 108854.

39. Adams, K.L.; Cronn, R.; Percifield, R.; Wendel, J.F. Genes duplicated by polyploidy show unequal contributions to the transcriptome and organ-specific reciprocal silencing. Proc. Natl. Acad. Sci. USA 2003, 100, 4649-4654.

40. Bowles, D.; Lim, E.K.; Poppenberger, B.; Vaistij, F.E. Glycosyltransferases of lipophilic small molecules. Annu. Rev. Plant Biol. 2006, 57, 567-597.

41. Ahrazem, O.; Rubio-Moraga, A.; Berman, J.; Capell, T.; Christou, P.; Zhu, C.; Gómez-Gómez, L. The carotenoid cleavage dioxygenase CCD2 catalyzing the synthesis of crocetin in spring crocuses and saffron is a plastidial enzyme. New Phytol. 2015, 209, 650-663.

42. Dar, A.A.; Choudhury, A.R.; Kancharla, P.K.; Arumugam, N. The FAD2 Gene in plants: Occurrence, regulation, and role. Front Plant Sci. 2017, 8, 1789.

43. Scotti-Campos, P.; Pais, I.P.; Ribeiro-Barros, A.I.; Martins, L.D.; Tomaz, M.A.; Rodrigues, W.P.; Campostrini, E.; Semedo, J.N.; Fortunato, A.S.; Martins, M.Q.; et al. Lipid profile adjustments may contribute to warming acclimation and to heat impact mitigation by elevated [CO2] in Coffea spp. Environ. Exp. Bot. 2019, 167, 103856.

44. Liu, D.; Gong, Q.; Ma, Y.; Li, P.; Li, J.; Yang, S.; Yuan, L.; Yu, Y.; Pan, D.; Xu, F.; et al. cpSecA, a thylakoid protein translocase subunit, is essential for photosynthetic development in Arabidopsis. J. Exp. Bot. 2010, 61, 1655-1669.

45. Kato, Y.; Kouso, T.; Sakamoto, W. Variegated tobacco leaves generated by chloroplast FtsH suppression: Implication of FtsH function in the maintenance of thylakoid membranes. Plant Cell Physiol. 2012, 53, 391-404.

46. Kaldenhoff, R.; Kai, L.; Uehlein, N. Aquaporins and membrane diffusion of CO2 in living organisms. Biochim. Biophys. Acta 2014, 1840, 1592-1595.

47. Karve, T.M.; Cheema, A.K. Small changes huge impact: The role of protein posttranslational modifications in cellular homeostasis and disease. J. Amino Acids 2011, 2011, 1-13.

48. Lyzenga, W.J.; Stone, S.L. Abiotic stress tolerance mediated by protein ubiquitination. J. Exp. Bot. 2012, 63, 599-616.

49. Erice, G.; Aranjuelo, I.; Irigoyen, J.J.; Sánchez-Díaz, M. Effect of elevated CO2, temperature and limited water supply on antioxidant status during regrowth of nodulated alfalfa. Physiol. Plant. 2007, 130, 33-45.

50. AbdElgawad, H.; Farfan-Vignolo, A.R.; Vos, D.; Asard, H. Elevated CO2 mitigates drought and temperature-induced oxidative stress differently in grasses and legumes. Plant Sci. 2015, 231, 1-10.

51. Watanabe, C.K.; Sato, S.; Yanagisawa, S.; Uesono, Y.; Terashima, I.; Noguchi, K. Effects of elevated CO2 on levels of primary metabolites and transcripts of genes encoding respiratory enzymes and their diurnal patterns in Arabidopsis thaliana: Possible relationships with respiratory rates. Plant Cell Physiol. 2014, 55, 341-357.

52. Vega, F.E.; Ziska, L.H.; Simpkins, A.; Infante, F.; Davis, A.P.; Rivera, J.A.; Barnaby, J.Y.; Wolf, J. Early growth phase and caffeine content response to recent and projected increases in atmospheric carbon dioxide in coffee (Coffea arabica and C. canephora). Sci Rep. 2020, 10, 5875.

53. Loreto, F.; Mannozzi, M.; Maris, C.; Nascetti, P.; Ferranti, F.; Pasqualini, S. Ozone quenching properties of isoprene and its antioxidant role in leaves. Plant Physiol. 2001, 126, 993-1000.

54. Sasaki, K.; Saito, T.; Lamsa, M.; Oksman-Caldentey, K.-M.; Suzuki, M.; Ohyama, K.; Muranaka, T.; Ohara, K.; Yazaki, K. Plants utilize isoprene emission as a thermo-tolerance mechanism. Plant Cell Physiol. 2007, 48, 1254-1262.

55. Belchínavarro, S.; Almagro, L.; Sabaterjara, A.B.; Fernández Pérez, F.; Bru, R.; Pedreño, M.A. Induction of trans-resveratrol and extracellular pathogenesis-related proteins in elicited suspension cultured cells of Vitis vinifera cv Monastrell. J. Plant Physiol. 2013, 170, 258-264.

56. Zhu, Q.G.; Gong, Z.Y.; Wang, M.M.; Li, X.; Grierson, D.; Yin, X.R.; Chen, K.S. A transcription factor network responsive to high CO2/hypoxia is in-volved in deastringency in persimmon fruit. J. Exp. Bot. 2018, 69, 2061-2070.

57. Wilkinson, M.J.; Monson, R.K.; Trahan, N.; Lee, S.; Brown, E.; Jackson, R.B.; Polley, H.W.; Fay, P.A.; Fall, R. Leaf isoprene emission rate as a function of atmospheric CO2 concentration. Glob. Chang. Biol. 2009, 15, 1189-1200.

58. Lee, S.B.; Suh, M.C. Advances in the understanding of cuticular waxes in Arabidopsis thaliana and crop species. Plant Cell Rep. 2015, 34, 557-572.

59. Gray, J.E.; Holroyd, G.H.; van der Lee, F.M.; Bahrami, A.R.; Sijmons, P.C.; Woodward, F.I.; Schuch, W.; Hetherington, A.M. The HIC signalling pathway links CO2 perception to stomatal development. Nature 2000, 408, 713.

60. Marsh, E.N.; Waugh, M.W. Aldehyde Decarbonylases: Enigmatic enzymes of hydrocarbon biosynthesis. ACS Catal. 2013, 3, 2515-2521.

61. Licausi, F.; Ohme-Takagi, M.; Perata, P. APETALA2/Ethylene responsive factor (AP2/ERF) transcription factors: Mediators of stress responses and developmental programs. New Phytol. 2013, 199, 639-649.

62. Müller, M.; Munné-Bosch, S. Ethylene response factors: A key regulatory hub in hormone and stress signalling. Plant Physiol. 2015, 169, 32-41.

63. Le Gall, H.; Philippe, F.; Domon, J.M.; Gillet, F.; Pelloux, J.; Rayon, C. Cell wall metabolism in response to abiotic stress. Plants 2015, 4, 112-166.

64. Uchida, A.; Jagendorf, A.T.; Hibino, T.; Takabe, T.; Takabe, T. Effects of hydrogen peroxide and nitric oxide on both salt and heat stress tolerance in rice. Plant Sci. 2002, 163, 515-523.

65. Ishibashi, Y.; Yamaguchi, H.; Yuasa, T.; Inwaya-Inoue, M.; Arima, S.; Zheng, S. Hydrogen peroxide spraying alleviates drought stress in soybean plants. J. Plant Physiol. 2011, 168, 1562-1567.

66. Hossain, M.A.; Fujita, M. Hydrogen peroxide priming stimulates drought tolerance in mustard (Brassica juncea L.) seedlings. Plant Gene Trait 2013, 4, 109-123.

67. Verbancic, J.; Lunn, J.E.; Stitt, M.; Persson, S. Carbon supply and the regulation of cell wall synthesis. Mol. Plant 2018, 11, 75-94.

68. Ezquer, I.; Salameh, I.; Colombo, L.; Kalaitzis, P. Plant cell walls tackling climate change: Biotechnological strategies to improve crop adaptations and photosynthesis in response to global warming. Plants 2020, 9, 212.

69. Ainsworth, E.A.; Long, S.P. What have we learned from 15 years of free air-CO2 enrichment (FACE)? A meta-analytic review of the responses of photosynthesis, canopy properties and plant production to rising CO2. New Phytol. 2005, 165, 351-372.

70. Yuan, H.; Zhang, J.; Nageswaran, D.; Li, L. Carotenoid metabolism and regulation in horticultural crops. Hortic. Res. 2015, 2, 15036.

71. Mondego, J.M.C.; Vidal, R.O.; Carazzolle, M.F. An EST-based analysis identifies new genes and reveals distinctive gene expression features of Coffea arabica and Coffea canephora. BMC Plant Biol. 2011, 11, 30.

72. Yuyama, P.M.; Júnior, O.R.; Ivamoto, S.T.; Silva Domingues, D.; Carazzolle, M.F.; Guimarães Pereira, G.A.; Charmetant, P.; Leroy, T.; Protasio Pereira, L.F. Transcriptome analysis in Coffea eugenioides, an Arabica coffee ancestor, reveals differentially expressed genes in leaves and fruits. Mol. Genet. Genomics 2016, 291, 323-336.

73. DaMatta, F.M.; Ramalho, J.D.C. Impacts of drought and temperature stress on coffee physiology and production: A review. Braz. J. Plant Physiol. 2006, 18, 55-81.

74. Ramalho, J.C.; DaMatta, F.M.; Rodrigues, A.P.; Scotti-Campos, P.; Pais, I.; Batista-Santos, P.; Partelli, F.L.; Ribeiro, A.; Lidon, F.C.; Leitão, A.E. Cold impact and acclimation response of Coffea spp. plants. Theor. Exp. Plant Physiol. 2014, 26, 5-18.

75. Colebrook, E.H.; Thomas, S.G.; Phillips, A.L.; Hedden, P. The role of gibberellin signalling in plant responses to abiotic stress. J. Exp. Biol. 2014, 217, 67-75.

76. Alves, G.; Torres, L.F.; Déchamp, E.; Breitler, J.C.; Joët, T.; Gatineau, F.; Andrade, A.C.; Bertrand, B.; Marraccini, P.; Etienne, H. Differential fine-tuning of gene expression regulation in coffee leaves by CcDREB1D promoter haplotypes under water deficit. J. Exp. Bot. 2017, 68, 3017-3031.

77. Torres, L.F.; Reichel, T.; Déchamp, E.; de Aquino, S.O.; Duarte, K.E.; Alves, G.S.C.; Silva, A.T.; Cotta, M.G.; Costa, T.S.; Diniz, L.E.C. Expression of DREB-Like genes in Coffea canephora and C. arabica subjected to various types of abiotic stress. Trop. Plant Biol. 2019, 12, 98-116.

78. Sakr, S.; Wang, M.; Dédaldéchamp, F.; Perez-Garcia, M.D.; Ogé, L.; Hamama, L.; Atanassova, R. The sugar-signaling hub: Overview of regulators and interaction with the hormonal and metabolic network. Int. J. Mol. Sci. 2018, 19, 2506.

79. Fukayama, H.; Sugino, M.; Fukuda, T.; Masumoto, C.; Taniguchi, Y.; Okada, M.; Sameshima, R.; Hatanaka, T.; Misoo, S.; Hasegawa, T.; et al. Gene expression profiling of rice grown in free air CO2 enrichment (FACE) and elevated soil temperature. Field Crop. Res. 2011, 121, 195-199.

80. Ghasemzadeh, A.; Jaafar, H.Z.E.; Rahmat, A. Elevated carbon dioxide increases contents of flavonoids and phenolic compounds, and antioxidant activities in Malaysian young ginger (Zingiber officinale Roscoe.) varieties. Molecules 2010, 15, 7907-7922.

81. Xu, Z.Z.; Shimizu, H.; Ito, S.; Yagasaki, Y.; Zou, C.J.; Zhou, G.S.; Zheng, Y. Effects of elevated CO2, warming and precipitation change on plant growth, photosynthesis and peroxidation in dominant species from North China grassland. Planta 2014, 239, 421-435.

82. Zinta, G.; AbdElgawad, H.; Domagalska, M.A.; Vergauwen, L.; Knapen, D.; Nijs, I.; Janssens, I.A.; Beemster, G.T.S.; Asard, H. Physiological, biochemical, and genome-wide transcriptional analysis reveals that elevated CO2 mitigates the impact of combined heat wave and drought stress in Arabidopsis thaliana at multiple organizational levels. Glob. Chang. Biol. 2014, 20, 3670-3685.

83. Pritchard, S.G.; Ju, Z.; Santen, E.; Qiu, J.; Weaver, D.B.; Prior, S.A.; Rogers, H.H. The influence of elevated CO2 on the activities of antioxidative enzymes in two soybean genotypes. Aust. J. Plant Physiol. 2000, 27, 1061-1068.

84. Ramalho, J.C.; Campos, P.S.; Teixeira, M.; Nunes, M.A. Nitrogen dependent changes in antioxidant systems and in fatty acid composition of chloroplast membranes from Coffea arabica L. plants submitted to high irradiance. Plant Sci. 1998, 135, 115-124.

85. Semedo, J.N.; Rodrigues, A.P.; Lidon, F.C.; Pais, I.P.; Marques, I.; Gouveia, D.; Armengaud, J.; Martins, S.; Semedo, M.C.; Silva, M.J.; et al. Intrinsic non-stomatal resilience to drought of the photosynthetic apparatus in Coffea spp. can be strengthened by elevated Air CO2. Tree Physiol. 2020, tpaa158.

86. Ijaq, J.; Chandrasekharan, M.; Poddar, R.; Bethi, N.; Sundararajan, V.S. Annotation and curation of uncharacterized proteins-challenges. Front. Genet. 2015, 6, 119.

87. Naveed, M.; Chaudhry, Z.; Ali, Z.; Amjad, M.; Zulfiqar, F.; Numan, A. Annotation and curation of hypothetical proteins: Prioritizing targets for experimental study. Adv. Life Sci. 2018, 5, 73-87.

88. Andrews, S. FastQC: A Quality Control Tool for High Throughput Sequence Data. 2010. Available online: http://www.bioinformatics.babraham.ac.uk/projects/fastqc (accessed on 12 November 2019).

89. Wingett, S.W.; Andrews, S. FastQ Screen: A tool for multi-genome mapping and quality control. F1000 Res. 2018, 7, 1338.

90. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114-2120.

91. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15-21.

92. Anders, S.; Pyl, P.T.; Huber, W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics 2015, 31, 166-169.

93. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R. 1000 genome project data processing subgroup. The sequence alignment/map (SAM) format and SAMtools. Bioinformatics 2009, 25, 2078-2079.

94. Pertea, G. Gffread: GFF/GTF Utility Providing Format Conversions, Region Filtering, FASTA Sequence Extraction and More. 2015. Available online: https://github.com/gpertea/gffread (accessed on 24 November 2019).

95. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018; Available online: https://www.R-project.org (accessed on 14 December 2019).

96. Trapnell, C.; Williams, B.A.; Pertea, G.; Mortazavi, A.; Kwan, G.; van Baren, M.J. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 2010, 28, 511-515.

97. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550.

98. Benjamini, Y.; Hochberg, Y. On the adaptive control of the false discovery rate in multiple testing with independent statistics. J. Educ. Behav. Sci. 2000, 25, 60-83.

99. Chen, H.; Boutros, P.C. VennDiagram: A package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinform. 2011, 12, 35.

100. Götz, S. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008, 36, 3420-3435.

101. Wickham, H. ggplot2: Elegant Graphics for Data Analysis; Springer: New York, NY, USA, 2016; ISBN 978-3-319-24277-4. Available online: https://ggplot2.tidyverse.org (accessed on 21 June 2020).

102. Binns, D.; Dimmer, E.; Huntley, R.; Barrell, D.; O'Donovan, C.; Apweiler, R. QuickGO: A web-based tool for Gene Ontology searching. Bioinformatics 2009, 25, 3045-3046.