A Transcriptomic Approach to Understanding the Combined Impacts of Supra-Optimal Temperatures and CO2 Revealed Different Responses in the Polyploid Coffea arabica and Its Diploid Progenitor C. canephora

International Journal of Molecular Sciences Article A Transcriptomic Approach to Understanding the Combined Impacts of Supra-Optimal Temperatures and CO2Revealed Different Responses in the Polyploid Coffea arabica and Its Diploid Progenitor C. canephora Isabel Marques 1,2,* , Isabel Fernandes 2, Octávio S. Paulo 2, Fernando C. Lidon 3, Fábio M. DaMatta 4, JoséC. Ramalho 1,3,* and Ana I. Ribeiro-Barros 1,3,*   Citation: Marques, I.; Fernandes, I.; Paulo, O.S.; Lidon, F.C.; DaMatta, F.M.; Ramalho, J.C.; Ribeiro-Barros, A.I. A Transcriptomic Approach to Understanding the Combined Impacts of Supra-Optimal Temperatures and CO2Revealed Different Responses in the Polyploid Coffea arabica and Its Diploid Progenitor C. canephora.Int. J. Mol. Sci. 2021,22, 3125. https:// doi.org/10.3390/ijms22063125 Academic Editor: Massimo Maffei Received: 1 February 2021 Accepted: 8 March 2021 Published: 18 March 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Plant-Environment Interactions and Biodiversity Lab (PlantStress&Biodiversity), Forest Research Centre (CEF), Instituto Superior de Agronomia (ISA), Universidade de Lisboa, 2784-505 Oeiras and Tapada da Ajuda, 1349-017 Lisboa, Portugal 2 Computational Biology and Population Genomics Group, Centre for Ecology, Evolution and Environmental Changes (cE3c), Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal; [email protected] (I.F.); [email protected] (O.S.P.) 3GeoBioSciences, GeoTechnologies and GeoEngineering (GeoBioTec), Faculdade de Ciências e Tecnologia (FCT), Universidade NOVA de Lisboa (UNL), 2829-516 Monte de Caparica, Portugal; [email protected] 4Departamento de Biologia Vegetal, Universidade Federal Viçosa (UFV), Viçosa 36570-090, MG, Brazil; [email protected] *Correspondence: [email protected] (I.M.); [email protected] (J.C.R.); [email protected] (A.I.R.-B.) Abstract: Understanding the effect of extreme temperatures and elevated air (CO 2 ) is crucial for mitigating the impacts of the coffee industry. In this work, leaf transcriptomic changes were evaluated in the diploid C. canephora and its polyploid C. arabica, grown at 25 ◦ C and at two supra-optimal temperatures (37 ◦ C, 42 ◦ C), under ambient (aCO 2 ) or elevated air CO 2 (eCO 2 ). Both species expressed fewer genes as temperature rose, although a high number of differentially expressed genes (DEGs) were observed, especially at 42 ◦ C. An enrichment analysis revealed that the two species reacted differently to the high temperatures but with an overall up-regulation of the photosynthetic machinery until 37 ◦ C. Although eCO 2 helped to release stress, 42 ◦ C had a severe impact on both species. A total of 667 photosynthetic and biochemical related-DEGs were altered with high temperatures and eCO 2 , which may be used as key probe genes in future studies. This was mostly felt in C. arabica, where genes related to ribulose-bisphosphate carboxylase (RuBisCO) activity, chlorophyll a-b binding, and the reaction centres of photosystems I and II were down-regulated, especially under 42 ◦ C, regardless of CO 2 . Transcriptomic changes showed that both species were strongly affected by the highest temperature, although they can endure higher temperatures (37 ◦ C) than previously assumed. Keywords: climate changes; coffee; elevated air (CO 2 ); functional analysis; high temperatures; leaf RNAseq; polyploidy; warming 1. Introduction Temperature and carbon dioxide (CO 2 ) are major drivers of climate change affecting global crop production [ 1 – 3 ]. It is therefore not surprising that many crop improvement programs focus on the development of climate-smart crops resilient to climate change [ 4 ]. Yet, breeding for tolerance to a single stress may be risky because in nature plants respond to simultaneous stresses, and increasing tolerance to one stress may be at the expense of resilience to another [ 5 , 6 ]. For instance, high temperatures cause physiological, biochemical, and molecular changes affecting key biological processes, such as photosynthesis, by reducing electron transport, NADPH and ATP synthesis, and ribulose-bisphosphate carboxylase Int. J. Mol. Sci. 2021,22, 3125. https://doi.org/10.3390/ijms22063125 https://www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2021,22, 3125 2 of 20 (RuBisCO) activity while increasing the production of H 2 O 2 [ 7 – 9 ]. At the other extreme, elevated CO 2 (eCO 2 ) was found to improve the physiological status and to mitigate the adverse effects of high temperatures in different species, increasing photosynthesis in most C3 plants [ 10 – 14 ]. Therefore, it is of utmost importance to identify the response pathways linked to several environmental changes, to best optimize crop improvement under the future estimated climate conditions. Coffee is among the most important agricultural products [ 15 ]. It is produced in about 80 tropical countries, with an annual production of around nine million tons of green beans and is estimated to involve 125 million people in its entire chain of value, with many smallholder farmers whose livelihoods are supported by this crop [ 16 – 18 ]. The Coffea genus comprises at least 125 species, although only two dominate the coffee trade: the allotetraploid C. arabica L. (Arabica coffee; 2n = 4x = 44) and one of its diploid ancestors, C. canephora Pierre ex A. Froehner (Robusta coffee; 2n = 2x = 22), with the former contributing around 60–65% of world production [ 19 ]. Coffea arabica is thought to be originated from a single natural polyploidization event between C. canephora and C. eugenioides that occurred in very recent evolutionary times (<50,000 years ago [ 20 ]). The quite low levels of genetic variation found in C. arabica when compared to its diploid progenitors are a concern in the context of environmental changes regarding the sustainability of this crop [ 20 ]. The predicted increase in global temperature could be catastrophic to C. arabica yields and quality, since some works consider this species more sensitive to elevated temperatures than C. canephora [ 15 ], with impacts being already felt on field plantations [ 21 , 22 ]. In fact, classical studies suggested that coffee photosynthesis would be particularly sensitive to temperatures above 20–25 ◦ C (e.g., [ 23 – 25 ]). It was accepted that temperature averages above 23 ◦ C accelerated the development and ripening of C. arabica fruits, often leading to the loss of coffee quality [ 26 ]. Continuous exposure to temperatures as high as 30 ◦C may restrain growth and can induce abnormalities such as the yellowing of leaves and growth of tumors at the base of the stem [ 27 ]. A relatively high temperature during blossoming, especially if associated with a prolonged dry season, may cause the abortion of flowers [ 26 ]. Nevertheless, selected cultivars under intensive management conditions have allowed C. arabica plantations to be spread to marginal regions with average temperatures as high as 24–25 ◦ C, with satisfactory yields, such as in the northeast of Brazil [ 28 , 29 ]. Furthermore, some experiments reported that if temperature rises gradually (from 24 ◦ C up to 33–35 ◦ C) and considering a sufficiently long acclimation time, C. arabica plants can increase their photosynthesis up to temperatures around 30 ◦ C, displaying the same efficiency as observed at 24 ◦ C [ 28 ]. Also, some genotypes can maintain full photosynthetic functioning at temperatures as high as 35 ◦ C [ 15 ] or even up to 39 ◦ C for a short number of days [ 14 , 30 ]. The maintenance of high photosynthetic efficiency suggests the absence of negative effects on the photosynthetic structures, which agrees to some extent with the absence of membrane permeability changes after a short-time exposure up to 50 ◦ C, as found in C. arabica [ 31 ]. For C. canephora, the optimum annual mean temperature ranges from 22 to 26 ◦ C [ 32 ], or from 24 to 30 ◦ C [ 33 ], but high temperatures can still be harmful, especially if the air is dry [ 34 ]. In this context, several breeding initiatives are being developed to adapt coffee crops to these environmental changes, looking back to genotypes that can be used to develop new resistant varieties while maintaining the highest cup quality ( https://worldcoffeeresearch.org/, accessed on 12 January 2021). In the last 20 years, a considerable amount of research has been devoted to environmental coffee physiology focusing on water relations and drought tolerance mechanisms, but little is still known on the mechanisms of tolerance to unfavorable temperatures, namely transcriptional and metabolic differences in response to warmer temperatures (see [ 35 ]). For instance, it has been observed that in Coffea, eCO 2 promotes a high plant vigor [ 13 ] and even crop yield [ 36 ], while improving tolerance to drought [ 37 – 39 ] and supra-optimal temperatures [ 14 , 37 , 40 , 41 ], with positive impacts on the physical and chemical traits of the coffee beans under supra-optimal temperatures, contributing to preserve its quality [ 42 ]. We previously showed that eCO 2 causes significant changes in the transcriptomic responses Int. J. Mol. Sci. 2021,22, 3125 3 of 20 of C. arabica and C. canephora, with differentially expressed genes (DEGs) being much more abundant in the diploid than in C. arabica [ 43 ]. Functional analysis also revealed that under eCO 2 ,C. canephora expressed DEGs associated to general biological processes and to a lower extent to abiotic stress responses while C. arabica showed an upregulation of genes linked to plant tolerance, oxidative stress control, aquaporins, reorganization of the lipid matrix membrane, chloroplast/thylakoid organization, and PS II repair. This suggested a protective role of eCO2under stressful conditions. In the present study we explored the transcriptome of C. arabica (cv. Icatu) and C. canephora (cv. Conilon Clone 153) leaves for a thorough understanding of how coffee plants regulate transcriptomic responses to elevated temperatures, and how eCO 2 can mitigate some deleterious effects at the physiological and metabolic level. For this, we have used a comparative approach using two coffee genotypes grown under a control temperature of 25/20 ◦ C (day/night) and then gradually exposed to increased temperatures up to 42/30 ◦ C in combination with ambient CO 2 (aCO 2 : 380 µ L L −1 ) or elevated CO 2 (eCO 2 : 700 µ L L −1 ) air. This analysis is complementary to (and advances) previous studies, in which the plants from these two genotypes showed a temperature tolerance up to 37/28 ◦C [ 30 ], together with an increased photosynthetic potential under eCO 2 [ 43 ], which might contribute to the relief of some negative effects imposed by harsh elevated temperatures. 2. Results 2.1. Overall Transcriptome Profiling and Mapping Statistics Stringent quality assessment and data filtering generated an average of 26.3 million clean reads, from an average 28.3 million raw reads. Overall, a high proportion of reads were aligned to the reference genome since only an average of 3.17% of reads were not mapped (Table S1). Statistical details for each replicate are depicted in Table S1. 2.2. General Patterns of Gene Expression The number of expressed genes varied widely between the control temperature at 25/20 ◦ C (25 ◦ C) and the two different supra-optimal temperatures of 37/28 ◦ C (37 ◦ C) and 42/30 ◦ C (42 ◦ C), either under 380 µ L L −1 aCO 2 or 700 µ L L −1 eCO 2 . Overall, fewer genes were expressed as temperatures increased, especially in combination with eCO 2 where the lowest number of expressed genes was observed (Figure 1A; Table S2). The principal coordinate analysis (PCoA) revealed a stronger transcriptomic response in the two species as a result of the highest temperature (Figure 1B). In fact, PC2 separated all plants under 42 ◦ C in the upper half of the axis from plants under 25 ◦ C and 37 ◦ C that were grouped below this axis (Figure 1B), irrespective of the CO 2 condition and the species involved. PC1 was able to separate and group all Conilon Clone 153 (CL153) plants under 25 ◦ C and 37 ◦ C in the lower right part of this axis, whereas Icatu plants showed a wider variation. However, it seems remarkable that at the highest temperature, the CO 2 condition become determinant with PC1 separating the plants under aCO 2 (left upper quadrant) from those under eCO2(right upper quadrant). 2.3. Differential Gene Expression Changes in Response to Supra-Optimal Temperatures An average of 11,170 and 11,833 DEGs were identified, respectively by DESeq2 and edgeR, with 93.34% identified by both tools. Detailed results of each method are depicted in Table S2. The number of DEGs varied from 9545 (37 ◦ C; eCO 2 ) to 13,134 (42 ◦ C; eCO 2 ) in Icatu and from 8240 (37 ◦ C; aCO 2 ) to 12,115 (42 ◦ C; eCO 2 ) in CL153, being consistently higher at the highest temperature in both species, especially when combined with eCO 2 (Figure 2). Under the same CO 2 , the majority of DEGs were shared between the two supra-optimal temperatures (Figure 2). Int. J. Mol. Sci. 2021,22, 3125 4 of 20 Figure 1. Differences in the patterns of gene expression. ( A ) Number of expressed genes in plants of C. arabica cv. Icatu and C. canephora cv. Conilon Clone 153 (CL153), grown either in ambient CO 2 (aCO 2 , 380 µ L L −1 ) or elevated CO 2 (eCO 2 , 700 µ L L −1 ) air at control temperature conditions (25/20 ◦ C, day/night; 25 ◦ C), a moderate supra-optimal temperature (37/28 ◦ C; 37 ◦ C), and an extreme supra-optimal temperature (42/30 ◦ C; 42 ◦ C). ( B ) Principal coordinate analysis (PCoA) of rlog transformed gene expression data generated by RNA-sequencing. Each treatment contains three biological replicates and is indicated by the colors depicted in A. The percentage of variance is indicated in each axis. Square symbols indicate Icatu, while circles indicate CL153. Detailed information is additionally given in Table S1. In Icatu, 14.6% (eCO 2 ) to 19.9% (aCO 2 ) of DEGs were specific of 37 ◦ C, while 30.0% (aCO 2 ) to 37.9% (eCO 2 ) were specific of 42 ◦ C (Figure 2). As such, a higher number of treatment-specific DEGs were reported with the rise of temperature, especially under eCO 2 where the highest number of DEGs were found (4980). This reflects a very specific response to the highest temperature that was further amplified by eCO2(Figure 2A,B). In comparison, CL153 showed a similar global % of specific DEGs with 14.4% to 16.7% of DEGs being expressed at 37 ◦ C, and 30.8% to 35.0% DEGs at 42 ◦ C (Figure 2), in eCO 2 and aCO 2 plants, respectively. However, although the absolute value of treatment-specific DEGs in CL153 was similar under aCO 2 (3703) and eCO 2 (3735), as temperatures rose its percentage was lower under eCO2than under aCO2(31% vs. 35%) (Figure 2C,D). Int. J. Mol. Sci. 2021,22, 3125 5 of 20 Figure 2. Transcriptional patterns among differentially expressed genes (DEGs) between the two supra-optimal temperatures of 37/28 ◦ C (37 ◦ C) and 42/30 ◦ C (42 ◦ C). Numbers represent the DEGs shared between treatments (blue), specific for 37 ◦ C (yellow) and specific for 42 ◦ C (orange), found in plants of Icatu (A,B) and CL153 (C,D) grown under 380 µL L−1aCO2or 700 µL L−1eCO2. The number of upand down-regulated DEGs in plants at 37 ◦ C and 42 ◦ C was nearly the same in both genotypes, except at 42 ◦ C in eCO 2 where DEGs were predominantly up-regulated in Icatu (Figure 3). The expression profiles of specific DEGs were particularly influenced by the highest temperature, especially in Icatu which displayed a higher degree of variation (Figure 4). By contrast, the expression profile of the two genotypes exhibited a small differentiation at 37 ◦C. Figure 3. The effect of the supra-optimal temperatures 37/28 ◦ C (37 ◦ C) and 42/30 ◦ C (42 ◦ C) on the number of upand down-regulated treatment-specific DEGs in plants of Icatu and CL153, grown in either 380 µL L−1aCO2or 700 µL L−1eCO2. Int. J. Mol. Sci. 2021,22, 3125 6 of 20 Figure 4. Clustered heat maps and dendrograms of the normalized log 2 fold change visualizing the expression of treatment-specific differentially expressed genes (DEGs) in Icatu and CL153 as a response to 37/28 ◦ C (37 ◦ C) and 42/30 ◦ C (42 ◦ C) temperatures under 380 µ L L −1 aCO 2 or 700 µL L−1eCO2 . Significant DEGs were filtered by log 2 fold change (FC) > |2| and the plotted values scaled by row using Z-scores. Hot colors represent up-regulated DEGs and cold colors represent down-regulated DEGs. Column color labels group comparisons by temperature treatments (yellow: 37 ◦C; orange: 42 ◦C). 2.4. Significantly Enriched GO Terms of Responsive DEGs in the Two Genotypes In both species, an average of 74% of DEGs were annotated with gene ontology (GO) terms using the functional annotation of the reference genome of C. canephora (Table S3). Temperature increase had a negative effect on enriched GO terms, especially under eCO2. Up-regulated DEGs at 37 ◦ C showed more enriched GO terms than at 42 ◦ C, under both CO 2 conditions (Figure 5A; Table S4). Within the same temperature treatment, eCO 2 triggered more enriched GO terms than aCO 2 . The most enriched GO term (almost 300 genes) was related to molecular functions (RNA binding, GO:0003723) and was found only at 42 ◦ C in eCO 2 . The remaining were less enriched categories with a predominance of GO terms related to photosynthesis (GO:0015979), thylakoid (GO:0009579), photosystem (GO:0009521), photosynthesis, light reaction (GO:0019684), and chlorophyll binding (GO:0016168) found only at 37 ◦ C. Up-regulated DEGs in 42 ◦ C Icatu plants were also enriched for heat shock protein binding (GO:0031072) but only under aCO 2 . In fact, at 42 ◦C in eCO2there was an enrichment in the folding and binding of proteins. Int. J. Mol. Sci. 2021,22, 3125 7 of 20 Figure 5. Over-representation analysis of gene ontology (GO) terms performed with gProfiler against the functional annotation of the Coffea canephora genome. Significantly (false discovery rate (FDR) < 0.01) enriched gene ontology (GO) terms among up-regulated ( A ) and down-regulated ( B ) differentially expressed genes (DEGs) in Icatu and CL153 were ranked by increasing log 2 fold change (FC), considering the effect of supra-optimal temperatures at 37/28 ◦ C (37 ◦ C), and 42/30 ◦ C (42 ◦ C) in plants grown in 380 µ L L −1 aCO 2 or 700 µ L L −1 eCO 2 . The absence of a treatment in the figure indicates that no enriched GO terms were found. GO terms are grouped by the main categories: biological process (GO:BP), molecular function (GO:MF), and cellular component (GO:CC). Counts indicate the number of DEGs annotated with each GO term and dots are colored by temperature treatment. Down-regulated DEGs featured 37 enriched categories (17 in Icatu and 20 in CL153), although all with a very small number of gene counts (Figure 5B; Table S5). At 37 ◦ C, down-regulated DEGs showed almost no enriched categories under aCO 2 while under Int. J. Mol. Sci. 2021,22, 3125 8 of 20 eCO 2 , Icatu and especially CL153 showed an enrichment in general activities as for instance, oxidoreductase activity acting on paired donors with incorporation or reduction of molecular oxygen (GO:0016705). At 42 ◦ C in aCO 2 , down-regulated DEGs linked with microtubule motor activity (GO:0003777) and microtubule binding (GO:0008017) were highly enriched in Icatu while CL153 showed an enrichment in the molecular functions related to transporter activity. In this supra-optimal temperature, eCO 2 triggered less enriched categories: Icatu plants were mostly enriched in two biological processes, the secondary metabolic process (GO:0019748), and the lignin catabolic process (GO:0046274), while CL153 plants showed an enrichment in the molecular functions linked to calcium ion binding (GO:0005509), xyloglucan:xyloglucosyl transferase activity (GO:0016762), sulfotransferase activity (GO:0008146), and oxidoreductase activity, acting on diphenols and related substances as donors, oxygen as an acceptor (GO:0016682). 2.5. Effect of Supra-Optimal Temperatures and eCO2on Photosynthetic and Other Biochemical-Related Responsive DEGs A total of 667 responsive DEGs associated with photosynthesis and other important metabolic functions in coffee leaves (e.g., antioxidant mechanisms, lipids, and respiratory pathways) were found in the two genotypes: 180 directly related with the photosynthetic pathway (122 to photosynthesis, 49 to the chlorophyll metabolic process, 9 to RuBisCO activity), 73 to antioxidant activities, 123 to lipid metabolism, and 291 directly related with the respiratory pathway (170 to cellular respiration, 30 to malate dehydrogenase activity, and 91 to pyruvate kinase activity) (Table S6). Fold changes of these DEGs varied widely but the most extreme values were always found at 42 ◦C. Overall, a similar number of photosynthetic DEGs were found in the two genotypes as a response to supra-optimal temperatures although with contrasting levels of expression. These DEGs were mostly up-regulated at 37 ◦ C but down-regulated at 42 ◦ C in Icatu (Figure 6A), while the majority were up-regulated in CL153 (Figures 6B and 7). The same pattern was found for DEGs involved in the chlorophyll metabolic process (Figures 6and 7) . In comparison, DEGs linked to RuBisCO activity were mostly down-regulated in Icatu plants under 42 ◦ C independently of CO 2 , while in CL153 they were always up-regulated (Figures 6and 7). More than half of the photosynthesis DEGs were involved in binding activities (Table S6). Notably, under 42 ◦ C and independently of the CO 2 conditions, DEGs related to the reaction centres of photosystems (PSs) I and II were down-regulated in Icatu, contrary to CL153 where they were up-regulated. Under these extreme conditions, Icatu also repressed genes involved in chlorophyll a-b binding and most PsbQ and PsbP genes, contrary to CL153 where they were generally up-regulated (Table S6). The remaining categories were mostly down-regulated, especially at 42 ◦ C under eCO 2 although differences were still found between the two genotypes (Figure 6A,B). For instance, DEGs involved in antioxidant activities and lipid metabolism were always more down-regulated in Icatu except at 37 ◦ C in eCO 2 , while in CL153 those DEGs were always up-regulated at 37 ◦ C, independently of CO 2 (Figure 6A,B). Overall, in Icatu plants under both CO 2 conditions, the number of up-regulated DEGs involved in cellular respiration was higher in 37 ◦ C plants than in the 42 ◦ C plants. In CL153, a high number of upregulated DEGs were found, but only in 37 ◦ C plants under aCO 2 . DEGs involved in pyruvate kinase (PK) and malate dehydrogenase (MDH) activity (involved in glycolysis) were mostly down-regulated in all treatments, except for MDH in 37 ◦ C-CL153 plants under aCO2. Int. J. Mol. Sci. 2021,22, 3125 9 of 20 Figure 6. Changes in the proportion of the regulation (%) and in the numbers (indicated in each bar) of differentially expressed genes (DEGs) related to photosynthesis and biochemical processes in Icatu ( A ) and CL153 ( B ) plants, as a response to 37/28 ◦ C (37 ◦ C) and 42/30 ◦ C (42 ◦ C) and grown in 380 µL L−1 aCO 2 or 700 µ L L −1 eCO 2 . The searched gene ontology (GO) biological processes included: “photosynthesis”, “chlorophyll metabolic process”, “ribulose-bisphosphate carboxylase activity” (RuBisCO), “antioxidant activity”, “lipid metabolic process (LOX, FAD)”, “cellular respiration”, “malate dehydrogenase activity”, and “pyruvate kinase activity”, as well as its direct child significant terms. Figure 7. Heatmap and dendrograms of the normalized log 2 fold change (FC) of photosynthesis-related significant differentially expressed genes (DEGs) of Icatu and CL153 plants as a response to 37/28 ◦ C (37 ◦ C) and 42/30 ◦ C (42 ◦ C) and grown in 380 µ L L −1 aCO 2 or 700 µ L L −1 eCO 2 . DEGs presented here are annotated with the gene ontology (GO) terms “photosynthesis”, “chlorophyll metabolism process”, and “RuBisCO”, according to the functional annotation of the Coffea canephora genome. The plotted values were scaled by row for improved visualization. Hot colors represent up-regulated DEGs, and cold colors represent down-regulated DEGs. Column color labels group comparisons by temperature treatment, while row color labels group genes by GO annotation. Int. J. Mol. Sci. 2021,22, 3125 16 of 20 the different treatments. To prevent highly differentially expressed genes from clustering together without considering their expression pattern, the log2fold change was scaled by gene across treatments (row Z-score). Supplementary Materials: The following are available online at https://www.mdpi.com/1422-006 7/22/6/3125/s1. Author Contributions: Conceptualization, A.I.R.-B. and J.C.R.; data curation, I.M., I.F., O.S.P., and J.C.R.; formal analysis, I.M., I.F., and O.S.P.; funding acquisition, J.C.R., A.I.R.-B., F.C.L., and F.M.D.; investigation, J.C.R., A.I.R.-B.; F.C.L., and F.M.D.; methodology, J.C.R. and A.I.R.-B.; project administration, J.C.R. and A.I.R.-B.; supervision, I.M., O.S.P., A.I.R.-B., and J.C.R.; writing—original draft, I.M., I.F., O.S.P., F.C.L., F.M.D., A.I.R.-B., and J.C.R.; writing—review & editing, I.M., A.I.R.-B., J.C.R., F.M.D., O.S.P., and I.F. 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 the 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/ASPAGR/31257/2017, and the research units UIDB/00239/2020 (CEF), UIDP/04035/2020 (GeoBioTec). 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.D., are also greatly acknowledged. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Raw reads have been deposited in the NCBI Sequence Read Archive, BioProject accession PRJNA630692. 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. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. References 1. Taub, D.R.; Miller, B.; Allen, H. Effects of elevated CO 2 on the protein concentration of food crops: A meta-analysis. Glob. Chang. Biol. 2007,14, 565–575. [CrossRef] 2. IPCC. Working Group III. Climate Change 2014. Mitigation of Climate Change. In Contributions to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Edenhofer, O., Pichs-Madruga, R., Sokona, Y., Farahani, F., Kadner, S., Seyboth, K., Adler, A., Baum, I., Brunner, S., Eickemeier, P., et al., Eds.; Cambridge University Press: New York, NY, USA, 2014; p. 1435. 3. Simpson, B.M. Preparing Smallholder Farm Families to Adapt to Climate Change. Pocket Guide 2: Managing Crop Resources; Catholic Relief Services: Baltimore, MD, USA, 2017. 4. Wheeler, T.; Von Braun, J. Climate Change Impacts on Global Food Security. Nat. Clim. Chang. 2013 ,341, 479–485. [CrossRef] [PubMed] 5. Mittler, R. Abiotic stress, the field environment and stress combination. Trends Plant Sci. 2006,11, 15–19. [CrossRef] [PubMed] 6. Atkinson, N.J.; Urwin, P.E. The interaction of plant biotic and abiotic stresses: From genes to the field. J. Exp. Bot. 2012 ,63, 3523–3543. [CrossRef] [PubMed] 7. Kurek, I.; Chang, T.K.; Bertain, S.M.; Madrigal, A.; Liu, L.; Lassner, M.W.; Zhu, G. Enhanced Thermostability of Arabidopsis Rubisco Activase Improves Photosynthesis and Growth Rates under Moderate Heat Stress. Plant Cell 2007 ,19, 3230–3241. [CrossRef] 8. Sharkey, T.D.; Zhang, R. High Temperature Effects on Electron and Proton Circuits of Photosynthesis. J. Integr. Plant Biol. 2010 ,52, 712–722. [CrossRef] 9. Bita, C.E.; Gerats, T. Plant tolerance to high temperature in a changing environment: Scientific fundamentals and production of heat stress-tolerant crops. Front. Plant Sci. 2013,4, 273. [CrossRef] 10. Iii, E.W.H.; Heckathorn, S.A.; Joshi, P.; Wang, D.; Barua, D. Interactive Effects of Elevated CO2 and Growth Temperature on the Tolerance of Photosynthesis to Acute Heat Stress in C3 and C4 Species. J. Integr. Plant Biol. 2008,50, 1375–1387. [CrossRef] 11. Leakey, A.D.B.; Xu, F.; Gillespie, K.M.; McGrath, J.M.; Ainsworth, E.A.; Ort, D.R. Genomic basis for stimulated respiration by plants growing under elevated carbon dioxide. Proc. Natl. Acad. Sci. USA 2009,106, 3597–3602. [CrossRef] Int. J. Mol. Sci. 2021,22, 3125 17 of 20 12. Prasad, P.V.V.; Vu, J.C.V.; Boote, K.J.; Allen, L.H. Enhancement in leaf photosynthesis and upregulation of Rubisco in the C4 sorghum plant at elevated growth carbon dioxide and temperature occur at early stages of leaf ontogeny. Funct. Plant Biol. 2009 , 36, 761–769. [CrossRef] 13. Ramalho, J.C.; Rodrigues, A.P.; Semedo, J.N.; Pais, I.P.; Martins, L.D.; Simões-Costa, M.C.; Leitão, A.E.; Fortunato, A.S.; BatistaSantos, P.; Palos, I.M.; et al. Sustained Photosynthetic Performance of Coffea spp. under Long-Term Enhanced [CO 2 ]. PLOS ONE 2013,8, e82712. [CrossRef] 14. Rodrigues, W.P.; Martins, M.Q.; Fortunato, A.S.; Rodrigues, A.P.; Semedo, J.N.; Simões-Costa, M.C.; Pais, I.P.; Leitao, A.E.; Colwell, F.; Goulão, L.; et al. Long-term elevated air [CO 2 ] strengthens photosynthetic functioning and mitigates the impact of supra-optimal temperatures in tropical Coffea arabica and C. canephora species. Glob. Chang. Biol. 2015 ,22, 415–431. [CrossRef] 15. 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. [CrossRef] 16. Osorio, N. The Global Coffee Crisis: A Threat to Sustainable Development. In Proceedings of the World Summit on Sustainable Development, Johannesburg, South Africa, 30 August 2002. 17. 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. [CrossRef] 18. Theory and Practice of Climate Adaptation. In Climate Change Management; Springer International Publishing: New York, NY, USA, 2018; pp. 465–477. ISBN 978-3-319-72874-2. 19. International Coffee Organization. ICO Annual Review 2010; ICO: London, UK, 2010. 20. 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. [CrossRef] 21. Van Der Vossen, H.; Bertrand, B.; Charrier, A. Next generation variety development for sustainable production of arabica coffee (Coffea arabica L.): A review. Euphytica 2015,204, 243–256. [CrossRef] 22. Venancio, L.P.; Filgueiras, R.; Mantovani, E.C.; Amaral, C.H.D.; Da Cunha, F.F.; Silva, F.C.D.S.; Althoff, D.; Dos Santos, R.A.; Cavatte, P.C. Impact of drought associated with high temperatures on Coffea canephora plantations: A case study in Espírito Santo State, Brazil. Sci. Rep. 2020,10, 1–21. [CrossRef] 23. Nunes, M.A.; Bierhuizen, J.F.; Ploegman, C. Studies on Productivity of Coffee: I. Effect of Light, Temperature and CO2 Concentration on Photosynthesis of Coffee Arabica. Acta Bot. Neerl. 1968,17, 93–102. [CrossRef] 24. Cannell, M.G. Crop physiological aspects of coffee bean yield: A review. Kenya Coffee 1976,41, 245–253. 25. Kumar, D.; Tieszen, L.L. Some aspects of photosynthesis and related processes in Coffea arabica L. Kenya Coffee 1976 ,41, 309–315. 26. Camargo, A.P. O clima e a cafeicultura no Brasil. Inf. Agropec. 1985,11, 13–26. 27. Franco, C.M. Influence of Temperature on Growth of Coffee Plant; Bulletin, No. 16; IBEC Research Institute: New York, NY, USA, 1958. 28. DaMatta, F.M. Ecophysiological constraints on the production of shaded and unshaded coffee: A review. Field Crop. Res. 2004 ,86, 99–114. [CrossRef] 29. Marie, L.; Abdallah, C.; Campa, C.; Courtel, P.; Bordeaux, M.; Navarini, L.; Lonzarich, V.; Bosselmann, A.S.; Turreira-García, N.; Alpizar, E.; et al. G × E interactions on yield and quality in Coffea arabica: New F1 hybrids outperform American cultivars. Euphytica 2020,216, 1–17. [CrossRef] 30. Dubberstein, D.; Lidon, F.C.; Rodrigues, A.P.; Semedo, J.N.; Marques, I.; Rodrigues, W.P.; Gouveia, D.; Armengaud, J.; Semedo, M.C.; Martins, S.; et al. Resilient and Sensitive Key Points of the Photosynthetic Machinery of Coffea spp. to the Single and Superimposed Exposure to Severe Drought and Heat Stresses. Front. Plant Sci. 2020,11, 1049. [CrossRef] 31. Gascó, A.; Nardini, A.; Salleo, S. Resistance to water flow through leaves of Coffea arabica is dominated by extra-vascular tissues. Funct. Plant Biol. 2004,31, 1161–1168. [CrossRef] 32. Matiello, J.B. CaféConillon: Como Plantar, Tratar, Colher, Preparar e Vender; MM Produções Gráficas: Rio de Janeiro, Brazil, 1998. 33. Willson, K.C. Coffee, Cocoa ant Tea; CAB International: Wallingford, UK, 1999. 34. Coste, R. Coffee: The Plant and the Product; MacMillan Press: London, UK, 1992. 35. De Oliveira, R.R.; Ribeiro, T.H.C.; Cardon, C.H.; Fedenia, L.; Maia, V.A.; Barbosa, B.C.F.; Caldeira, C.F.; Klein, P.E.; Chalfun-Junior, A. Elevated Temperatures Impose Transcriptional Constraints and Elicit Intraspecific Differences Between Coffee Genotypes. Front. Plant Sci. 2020,11, 1113. [CrossRef] 36. DaMatta, F.M.; Rahn, E.; Läderach, P.; Ghini, R.; Ramalho, J.C. Why could the coffee crop endure climate change and global warming to a greater extent than previously estimated? Clim. Chang. 2018,152, 167–178. [CrossRef] 37. 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. [CrossRef] 38. Avila, R.T.; de Almeida, W.L.; Costa, L.C.; Machado, K.L.; Barbosa, M.L.; de Souza, R.P.; Martino, P.B.; Juárez, M.A.; Marçal, D.M.; Martins, S.C.; et al. Elevated air [CO2] improves photosynthetic performance and alters biomass accumulation and partitioning in drought-stressed coffee plants. Environ. Exp. Bot. 2020,177, 104137. [CrossRef] 39. Semedo, J.N.; Rodrigues, A.P.; Lidon, F.C.; Pais, I.P.; Marques, I.; Gouveia, D.; Armengaud, J.; Silva, M.J.; Martins, S.; Semedo, M.C.; et al. Intrinsic non-stomatal resilience to drought of the photosynthetic apparatus in Coffea spp. is strengthened by elevated air [CO2]. Tree Physiol. 2020,2021. [CrossRef] Int. J. Mol. Sci. 2021,22, 3125 18 of 20 40. Martins, L.D.; Tomaz, M.A.; Lidon, F.C.; DaMatta, F.M.; Ramalho, J.C. Combined effects of elevated [CO2] and high temperature on leaf mineral balance in Coffea spp. plants. Clim. Chang. 2014,126, 365–379. [CrossRef] 41. Martins, M.Q.; Rodrigues, W.P.; Fortunato, A.S.; Leitão, 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. [CrossRef] 42. Ramalho, J.C.; Pais, I.P.; Leitão, A.E.; Guerra, M.; Reboredo, F.H.; Máguas, C.M.; Carvalho, M.L.; Scotti-Campos, P.; Ribeiro-Barros, A.I.; Lidon, F.J.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. [CrossRef] 43. Marques, I.; Fernandes, I.; David, P.H.; Paulo, O.S.; Goulao, L.F.; Fortunato, A.S.; Lidon, F.C.; DaMatta, F.M.; Ramalho, J.C.; Ribeiro-Barros, A.I. Transcriptomic Leaf Profiling Reveals Differential Responses of the Two Most Traded Coffee Species to Elevated [CO2]. Int. J. Mol. Sci. 2020,21, 9211. [CrossRef] [PubMed] 44. Dusenge, M.E.; Duarte, A.G.; Way, D.A. Plant carbon metabolism and climate change: Elevated CO2 and temperature impacts on photosynthesis, photorespiration and respiration. New Phytol. 2018,221, 32–49. [CrossRef] [PubMed] 45. Vico, G.; Way, D.A.; Hurry, V.; Manzoni, S. Can leaf net photosynthesis acclimate to rising and more variable temperatures? Plant Cell Environ. 2019,42, 1913–1928. [CrossRef] 46. Prasch, C.M.; Sonnewald, U. Signaling events in plants: Stress factors in combination change the picture. Environ. Exp. Bot. 2015 , 114, 4–14. [CrossRef] 47. Bertrand, B.; Bardil, A.; Baraille, H.; Dussert, S.; Doulbeau, S.; Dubois, E.; Severac, D.; Dereeper, A.; Etienne, H. The Greater Phenotypic Homeostasis of the Allopolyploid Coffea arabica Improved the Transcriptional Homeostasis Over that of Both Diploid Parents. Plant Cell Physiol. 2015,56, 2035–2051. [CrossRef] 48. Zhang, X.; Rerksiri, W.; Liu, A.; Zhou, X.; Xiong, H.; Xiang, J.; Chen, X.; Xiong, X. Transcriptome profile reveals heat response mechanism at molecular and metabolic levels in rice flag leaf. Gene 2013,530, 185–192. [CrossRef] 49. Liu, Z.; Xin, M.; Qin, J.; Peng, H.; Ni, Z.; Yao, Y.; Sun, Q. Temporal transcriptome profiling reveals expression partitioning of homeologous genes contributing to heat and drought acclimation in wheat (Triticum Aestivum L.). BMC Plant Biol. 2015 ,15, 1–20. [CrossRef] 50. Fernandes, J.; Morrow, D.J.; Casati, P.; Walbot, V. Distinctive transcriptome responses to adverse environmental conditions inZea maysL. Plant Biotechnol. J. 2008,6, 782–798. [CrossRef] [PubMed] 51. Ma, H.; Liu, M. The microtubule cytoskeleton acts as a sensor for stress response signaling in plants. Mol. Biol. Rep. 2019 ,46, 5603–5608. [CrossRef] [PubMed] 52. Hasanuzzaman, M.; Nahar, K.; Alam, M.; Roychowdhury, R.; Fujita, M. Physiological, Biochemical, and Molecular Mechanisms of Heat Stress Tolerance in Plants. Int. J. Mol. Sci. 2013,14, 9643–9684. [CrossRef] [PubMed] 53. Stratilová, B.; Kozmon, S.; Stratilová, E.; Hrmova, M. Plant Xyloglucan Xyloglucosyl Transferases and the Cell Wall Structure: Subtle but Significant. Molecules 2020,25, 5619. [CrossRef] 54. Perruc, E.; Charpenteau, M.; Ramirez, B.C.; Jauneau, A.; Galaud, J.P.; Ranjeva, R.; Ranty, B. A novel calmodu-lin-binding protein functions as a negative regulator of osmotic stress tolerance in Arabidopsis thaliana seedlings. Plant J. 2014 ,38, 410–420. [CrossRef] 55. Hirschmann, F.; Krause, F.; Papenbrock, J. The multi-protein family of sulfotransferases in plants: Composition, occurrence, substrate specificity, and functions. Front. Plant Sci. 2014,5, 556. [CrossRef] 56. Agrawal, A.A.; Hastings, A.P.; Johnson, M.T.J.; Maron, J.L.; Salminen, J.-P. Insect Herbivores Drive Real-Time Ecological and Evolutionary Change in Plant Populations. Science 2012,338, 113–116. [CrossRef] 57. Züst, T.; Rasmann, S.; Agrawal, A.A. Growth-defense tradeoffs for two major anti-herbivore traits of the common milkweed Asclepias syriaca. Oikos 2015,124, 1404–1415. [CrossRef] 58. Xie, M.; Zhang, J.; Tschaplinski, T.J.; Tuskan, G.A.; Chen, J.G.; Muchero, W. Regulation of Lignin Biosynthesis and Its Role in Growth-Defense Tradeoffs. Front. Plant Sci. 2018,9, 1427. [CrossRef] 59. United Nations. Transforming our World: The 2030 Agenda for Sustainable Development. 2015. Available online: https: //www.un.org/ga/search/view_doc.asp?symbol=A/RES/70/1&Lang=E (accessed on 17 July 2020). 60. Pais, I.P.; Reboredo, F.H.; Ramalho, J.C.; Pessoa, M.F.; Lidon, F.C.; Silva, M.M. Potential Impacts of Climate Change on Agriculture: A Review. Emir. J. Food Agric. 2020,32, 397–407. [CrossRef] 61. Rivero, R.M.; Mestre, T.C.; Mittler, R.; Rubio, F.; Garcia-Sanchez, F.; Martinez, V. The combined effect of salinity and heat reveals a specific physiological, biochemical and molecular response in tomato plants. Plant Cell Environ. 2013,37, 1059–1073. [CrossRef] 62. Shaar-Moshe, L.; Blumwald, E.; Peleg, Z. Unique Physiological and Transcriptional Shifts under Combinations of Salinity, Drought, and Heat. Plant Physiol. 2017,174, 421–434. [CrossRef] 63. Sage, R.F.; Kubien, D.S. The temperature response of C3and C4photosynthesis. Plant Cell Environ. 2007 ,30, 1086–1106. [CrossRef] [PubMed] 64. Way, D.A.; Oren, R.; Kroner, Y. The space-time continuum: The effects of elevated CO 2 and temperature on trees and the importance of scaling. Plant Cell Environ. 2015,38, 991–1007. [CrossRef] 65. Boisvenue, C.; Running, S.W. Impacts of climate change on natural forest productivity—evidence since the middle of the 20th century. Glob. Chang. Biol. 2006,12, 862–882. [CrossRef] Int. J. Mol. Sci. 2021,22, 3125 19 of 20 66. Chovancek, E.; Zivcak, M.; Botyanszka, L.; Hauptvogel, P.; Yang, X.; Misheva, S.; Hussain, S.; Brestic, M. Transient Heat Waves May Affect the Photosynthetic Capacity of Susceptible Wheat Genotypes Due to Insufficient Photosystem I Photoprotection. Plants 2019,8, 282. [CrossRef] [PubMed] 67. Nickelsen, J.; Rengstl, B. Photosystem II Assembly: From Cyanobacteria to Plants. Annu. Rev. Plant Biol. 2013 ,64, 609–635. [CrossRef] 68. Järvi, S.; Suorsa, M.; Aro, E.M. Photosystem II repair in plant chloroplasts: Regulation, assisting proteins and shared components with photosystem II biogenesis. Biochim. Biophys. Acta Bioenerg. 2015,1847, 900–909. [CrossRef] [PubMed] 69. Combes, M.C.; Cenci, A.; Baraille, H.; Bertrand, B.; Lashermes, P. Homeologous Gene Expression in Response to Growing Temperature in a Recent Allopolyploid (Coffea Arabica L.). J. Hered. 2011,103, 36–46. [CrossRef] [PubMed] 70. Ramalho, J.C.; Campos, P.S.; Teixeira, M.; Nunes, M. Nitrogen dependent changes in antioxidant system and in fatty acid composition of chloroplast membranes from Coffea arabica L. plants submitted to high irradiance. Plant Sci. 1998 ,135, 115–124. [CrossRef] 71. Fortunato, A.S.; Lidon, F.C.; Batista-Santos, P.; Leitão, A.E.; Pais, I.P.; Ribeiro, A.I.; Ramalho, J.C. Biochemical and molecular characterization of the antioxidative system of Coffea sp. under cold conditions in genotypes with contrasting tolerance. J. Plant Physiol. 2010,167, 333–342. [CrossRef] [PubMed] 72. Sdiri, S.; Rambla, J.L.; Besada, C.; Granell, A.; Salvador, A. Changes in the volatile profile of citrus fruit submitted to postharvest degreening treatment. Postharvest Biol. Technol. 2017,133, 48–56. [CrossRef] 73. 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. [CrossRef] 74. Moat, J.; Gole, T.W.; Davis, A.P. Least concern to endangered: Applying climate change projections profoundly influences the extinction risk assessment for wild Arabica coffee. Glob. Chang. Biol. 2018,25, 390–403. [CrossRef] 75. Vieira, A.; Silva, D.N.; Várzea, V.; Paulo, O.S.; Batista, D. Genome-Wide Signatures of Selection in Colletotrichum kahawae Reveal Candidate Genes Potentially Involved in Pathogenicity and Aggressiveness. Front. Microbiol. 2019 ,10, 1374. [CrossRef] [PubMed] 76. Zullo, J., Jr.; Pinto, H.S.; Assad, E.D.; Ávila, A.M.H. Potential for growing Arabica coffee in the extreme south of Brazil in a warmer world. Clim. Chang. 2011,109, 535–548. [CrossRef] 77. 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. [CrossRef] 78. Davis, A.P.; Chadburn, H.; Moat, J.; O’Sullivan, R.; Hargreaves, S.; Nic Lughadha, E. High extinction risk for wild coffee species and implications for coffee sector sustainability. Sci. Adv. 2019,5, eaav3473. [CrossRef] [PubMed] 79. Ramalho, J.C.; Fortunato, A.S.; Goulão, L.; Lidon, F.C. Cold-induced changes in mineral content in leaves of Coffea spp. Identification of descriptors for tolerance assessment. Biol. Plant. 2013,57, 495–506. [CrossRef] 80. 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 15 June 2020). 81. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014 ,30, 2114–2120. [CrossRef] 82. Wingett, S.W.; Andrews, S. FastQ Screen: A tool for multi-genome mapping and quality control. F1000Research 2018 ,7, 1338. [CrossRef] 83. 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–1184. [CrossRef] 84. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013,29, 15–21. [CrossRef] [PubMed] 85. Anders, S.; Pyl, P.T.; Huber, W. HTSeq: A Python framework to work with high-throughput sequencing data. Bioinformatics 2015 , 31, 166–169. [CrossRef] [PubMed] 86. 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. [CrossRef] 87. Pertea, G. 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). 88. 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. [CrossRef] [PubMed] 89. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2009,26, 139–140. [CrossRef] 90. Benjamini, Y.; Hochberg, Y. On the Adaptive Control of the False Discovery Rate in Multiple Testing with Independent Statistics. J. Educ. Behav. Stat. 2000,25, 60–83. [CrossRef] 91. Hunter, J.D. Matplotlib: A 2D Graphics Environment. Comput. Sci. Eng. 2007,9, 90–95. [CrossRef] 92. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2018. Int. J. Mol. Sci. 2021,22, 3125 20 of 20 93. Raudvere, U.; Kolberg, L.; Kuzmin, I.; Arak, T.; Adler, P.; Peterson, H.; Vilo, J. g:Profiler: A web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019,47, W191–W198. [CrossRef] [PubMed] 94. Supek, F.; Bošnjak, M.; Škunca, N.; Šmuc, T. REVIGO Summarizes and Visualizes Long Lists of Gene Ontology Terms. PLOS ONE 2011,6, e21800. [CrossRef] [PubMed] 95. Wickham, H. Elegant Graphics for Data Analysis; Springer International Publishing: New York, NY, USA, 2016.