Multi-Approach Analysis Reveals Local Adaptation in a Widespread Forest Tree of Reunion Island.

Detecting processes of local adaptation in forest trees and identifying environmental selective drivers are of primary importance for forest management and conservation. Transplant experiments, functional genomics and population genomics are complementary tools to efficiently characterize heritable phenotypic traits and to decipher the genetic bases of adaptive traits. Using an integrative approach combining phenotypic assessment in common garden, transcriptomics and landscape genomics, we investigated leaf adaptive traits in Coffea mauritiana, a forest tree endemic to Reunion Island. Eight populations of C. mauritiana originating from sites with contrasted environmental conditions were sampled in common garden to assess several leaf morphological traits, to analyze the leaf transcriptome and leaf cuticular wax composition. The relative alkane content of cuticular waxes was significantly correlated with major climatic gradients, paving the way for further transcriptome-based analyses. The expression pattern of cuticle biosynthetic genes was consistent with a modulation of alkane accumulation across the population studied, supporting the hypothesis that the composition of cuticular wax is involved in the local adaptation of C. mauritiana. Association tests in landscape genomics performed using RNA-seq-derived single-nucleotide polymorphisms revealed that genes associated with cell wall remodeling also likely play an adaptive role. By combining these different approaches, this study efficiently identified local adaptation processes in a non-model species. Our results provide the first evidence for local adaptation in trees endemic to Reunion Island and highlight the importance of cuticle composition for the adaptation of trees to the high evaporative demand in warm climates. � The Author(s) 2020. Published by Oxford University Press on behalf of Japanese Society of Plant Physiologists.

Accession number: The nucleotide sequences reported in this article have been submitted to the European Nucleotide Archive under accession numbers SAMEA5423234 to SAMEA5423258 (study no PRJEB31715).

Introduction

Forest ecosystems are home of most of the world’s terrestrial biodiversity. However, most of forest tree species are currently facing new ecological constraints associated with global climate change. Identifying the processes involved in forest tree adaptation is therefore of prime importance to ensure appropriate forest management and to mitigate potential negative impacts of climate change (Aitken and Whitlock 2013, Jaramillo-Correa et al. 2015, Razgour et al. 2019). In particular, local adaptation arises when populations experience contrasted environmental conditions: selection is expected to favor traits that increase the fitness of individuals in their local environment. It has important implications for response to global change. For instance, when environmental conditions are changing, pre-adapted variants can invade the local population through migration and selection, allowing for a relatively quick evolutionary response (Aitken and Whitlock 2013). Local adaptation appears more and more to be a major factor to be taken into account in forest ecosystem rehabilitation and restoration programmes (Thomas et al. 2014).

Local adaptation in forest trees has traditionally been investigated through in situ characterization, provenance experiments (i.e. common garden experiments), reciprocal transplants and population genetics (Kawecki and Ebert 2004, Leinonen ). By growing accessions from different populations in similar environmental conditions, it is possible to study the genetic bases of adaptive traits while limiting the confounding effects of the environment from which the population originated (Gonzalez-Martinez et al. 2006, de Villemereuil et al. 2016). Putative adaptive traits can be identified by correlating phenotypic traits assessed in a common environment with the climatic conditions at the site of origin of the population concerned (Gonzalez-Martinez et al. 2006, de Villemereuil et al. 2016). However, while providing evidence for adaptive phenotypic traits, studying phenotypic polymorphism in a common environment does not provide information on the specific genes underlying adaptive traits (Gonzalez-Martinez et al. 2006, de Villemereuil et al. 2016). Advances in molecular biology have provided new approaches, such as functional genomics and population genomics, that can more accurately characterize the molecular processes underlying the local adaptation of trees (Gonzalez-Martinez et al. 2006, Sork et al. 2013). Despite its high potential (De Kort et al. 2014, de Villemereuil et al. 2016), few studies have used a combination of classical and modern approaches to find signatures of local adaptation.

Genome-wide molecular markers make it possible to screen genomes to identify differentiated regions (i.e. outlier loci) that are putatively under natural selection (Manel et al. 2010, Narum et al. 2013, De Kort et al. 2014). As single-nucleotide polymorphisms (SNPs) are abundant and spread over the genome, they are considered as well suited to study local adaptation. SNPs can be associated with non-synonymous substitutions in coding regions, changes in the transcription level of functional genes and in the regulation processes underlying adaptive phenotypes (Gonzalez-Martinez et al. 2006, Sork et al. 2013). Landscape genomics approaches aim to test for associations between putative adaptive loci (e.g. SNPs) and environmental descriptors of the species habitat while accounting for neutral patterns that affect allelic frequencies (i.e. genetic structure and demographic history) (Joost et al. 2007, Coop et al. 2010, Manel et al. 2010, Sork et al. 2013). Thus, landscape genomics approaches not only provide candidate loci for adaptation but also make it possible to identify the ecological pressures most likely responsible for local adaptation processes. Like most genomics approaches, landscape genomics approaches only provide partial information on the genetic basis of local adaptation. In fact, many adaptive traits are likely governed by multiple loci, each contributing only modestly to the adaptive trait concerned (Pritchard et al. 2010). As landscape genomics approaches often rely on identifying individual loci that look unusual, these approaches are not sufficiently powerful to identify adaptive processes governed by multiple loci (Pritchard et al. 2010). Consequently, approaches focusing on groups of genes that underpin biological processes, such as metabolic pathways, usefully complement the results of landscape genomics studies. Focusing on a well-documented candidate pathway is particularly useful in the case of species with limited genomic knowledge, as is the case of a number of tropical forest trees.

Among candidate pathways for local adaptation, changes in leaf cuticular composition have been reported as essential in adapting to a variety of environmental constraints (Lewandowska et al. 2020). Plant cuticle is composed of a cutin polyester matrix and intracuticular and epicuticular waxes that form a hydrophobic surface that protects the plants against different threats, including dehydration, ultraviolet radiation and pathogen attacks (Yeats and Rose 2013). Cuticular wax consists of derivatives of very-long-chain fatty acid (VLCFA) and its biosynthesis is accomplished in several major steps (Kunst and Samuels 2009, Bernard and Joub�s 2013). In plants, de novo fatty acid (FA) synthesis occurs in the plastid and consists of sequential reactions that add two-carbon units to the elongating acyl chain, which is conjugated to acyl carrier protein (ACP) during this cyclic process. Nascent C16 and C18 FAs are then transferred to the endoplasmic reticulum and activated into acyl-CoA by long-chain acyl-CoA synthetases (LACS). Fatty acyl-CoAs (C16 and C18) are then elongated to wax precursors (VLCFAs) with C26 to C34 chains by a repeating reaction process via the fatty acid elongase (FAE) complex. Following elongation, the conversion of VLCFAs to the different wax components involves two pathways: the decarbonylation pathway produces alkanes, secondary alcohols and ketones, while the acyl-reduction pathway produces primary alcohols and wax esters (Kunst and Samuels 2009, Bernard and Joub�s 2013, Borisjuk et al. 2014). Changes in the leaf cuticular wax load and composition under drought conditions have been reported in many species, including Arabidopsis (Kosma et al. 2009), crops, such as cotton, wheat and sesame (Bondada et al. 1996, Kim et al. 2007, Li et al. 2019), and trees (Cameron et al. 2006, Le Provost et al. 2013). In most cases, the most significant change is an increase in the level of very-long-chain alkanes under water stress (Bondada et al. 1996, Kim et al. 2007, Kosma et al. 2009, Li et al. 2019). How a water deficit affects the expression of genes involved in wax alkane biosynthesis has been largely elucidated in Arabidopsis (Kosma et al. 2009, Bourdenx et al. 2011, Kim et al. 2019). In forest trees, alkane composition and abundance appeared to vary with the environment, i.e. latitude, elevation, illumination and aridity (Hoffmann et al. 2013, Mitić et al. 2016, Feakins et al. 2016, Lykholat et al. 2018, Teunissen van Manen et al. 2019). Wax biosynthesis is well-known for its plasticity, but some studies suggest that it may also result from heritable factors linked with the geographic origin and the environmental conditions of forest trees (Mitić et al. 2016, Li et al. 2020).

Coffea mauritiana is an endemic tree species from the tropical rain forest of Reunion and Mauritius islands. Related to the cultivated coffee tree species (C. arabica and C. canephora), it belongs to the Rubiaceae family. Due to human activities, only small patches of relict forest remain in Mauritius, whereas most of this forest is well preserved in Reunion Island where C. mauritiana is widely distributed (270–1,500 m a.s.l.) (Dulloo et al. 1999), in the form of small patches of individuals. According to genetic studies, C. mauritiana individuals mainly cluster according to the site of plant sampling (source populations), with further grouping in genetic clusters showing regional distribution (Garot et al. 2019a). During the Last Glaciation Maximum, the size of C. mauritiana population was reduced (Garot et al. 2019b) and the distribution of this species was likely restricted to the rainy windward side of the island. The current species distribution appears to be associated with a recent expansion of the species across the island (Garot et al. 2019b) and includes both gradients of temperature (13–24�C in mean annual temperature) and rainfall (from 1,100 to up to 5,000 mm in median annual precipitation) (Supplementary Fig. S1). The availability of ecological opportunity in young insular environments and intrinsic and extrinsic factors promoting the acquisition or retention of genetic variation have long been recognized as important considerations in assessing the potential of plant lineages for island colonization and diversification (Baldwin 2019). Local adaptation may emerge when natural selection by environmental factors is stronger than the homogenizing effect of gene flow among populations (Izuno et al. 2017). A species, such as C. mauritiana, occupying a wide range of ecological niches in Reunion Island provides a means to study intraspecific genetic diversification, in which the evidence and trajectories of adaptive evolution are likely to be preserved.

By examining this system, i.e. C. mauritiana in Reunion Island, we aimed to understand how a tree species adapts to contrasted climatic conditions on the spatial scale of a small oceanic island. Specifically, we tested whether leaf characteristics contribute to the local adaptation of C. mauritiana and whether there are genomic signatures of differentiation underlying adaptive traits among populations. To address these questions, we used a combination of classic and modern approaches, involving phenotypic measurements of leaf traits, transcriptomic and genomic analyses. The experimental design included eight tree populations grown in common garden (collection) and trees of in situ populations located across Reunion Island and covering most of the environmental distribution of the species. Phenotypic traits, including leaf morphological traits and cuticular wax composition (for common garden trees only), were measured. First insights into local adaptation were obtained by correlating the phenotypic traits assessed in common garden experiments with climatic factors at the original population sites. An RNA-seq experiment was also performed on these populations. Focusing on a candidate pathway, the expression of cuticle biosynthetic genes was retrieved from RNA-seq experiment and correlated with the proportion of alkanes present in the leaf cuticular waxes. In addition, landscape genomics studies were performed using RNA-seq-derived SNP to detect candidate genes under selection that indicate possible adaptation to climatic gradients in this species.

Results

Assessing correlations among climatic parameters

Possible correlations between the climatic parameters of the C. mauritiana population collection sites were examined. Pearson’s correlation coefficients and principal component analysis (PCA) analysis revealed that climatic parameters related to temperature [annual normal daily average temperature (Tmean), annual normal maximum daily temperature (Tmax) and annual normal minimum daily temperature (Tmin)] formed a cluster of highly correlated parameters (|r| > 0.8) separated from the climatic parameters related to precipitation [annual normal precipitations (Pmean) and annual normal precipitations during the dry season (July–September (Pmean_JAS)] (Fig.�1, Supplementary Table S1). Regarding evaporative demand, the annual normal potential evapotranspiration (PET) was highly correlated with temperature parameters (|r| =0.80 with Tmin; |r| =0.89 with Tmean; |r| =0.92 with Tmax) but showed relatively low correlation with precipitation parameters (|r| =0.49 with Pmean_JAS and |r| =0.17 with Pmean). Finally, the annual normal temperature amplitude (Trange) showed relatively low correlation with all other climatic parameters.

Fig. 1

Principal component analysis of the seven climatic factors. Percentage of explained variance of each axis is given in parentheses.

Selection of phenotypic traits

Descriptors of leaf morphology were derived from length and mass measurements performed on C. mauritiana populations sampled in situ and taken from the collection. Among morphological traits evaluated in the collection, basic measurements of leaf length, width and mass, as well as leaf area (LA) were highly correlated (Fig.�2a). Since they showed the lowest correlations, specific leaf area (SLA), leaf dry matter content (LDMC), leaf thickness (LTE) and LA were retained as variables for subsequent analysis.

Fig. 2

PCA of the phenotypic traits assessed in the collection. The percentage of explained variance of each axis is given in parentheses. (a) Results of the PCA of the eight leaf morphological traits. (b) Results of PCA of the four biochemical traits.

Among the biochemical traits of the cuticular wax evaluated in the collection, the total proportions of primary alcohols (OH) and alkane (Alk) were negatively correlated (Fig.�2b). Furthermore, primary alcohols and alkanes were the main classes of cuticular wax found in C. mauritiana leaves, together accounting for 82% (�7%) of the total wax load, while VLCFAs represented only minor components. More specifically, C28 and C30 chain lengths predominated in the primary alcohols (Fig.�3a), while the main alkanes had C29 and C31 chain lengths (Fig.�3b). Among the trees from the different populations, the total proportion of primary alcohols was highly negatively correlated (P-value = 3.71 � 10−4; R� = 90%) with that of alkanes (Fig.�3c). To limit redundancy, only the total amount of cuticular wax (wax.load) and the total proportion of Alk and FA were retained for subsequent analyses.

Fig. 3

Cuticular wax profiles of C. mauritiana from 23 accessions sampled in the collection. (a) Proportion of alcohols in total alcohol content. Stars represent outlier values. (b) Proportion of alkanes in total alkane content. Stars represent outlier values. (c) Linear regression among populations between the proportions of primary alcohols and alkanes.

Exploring phenotypic polymorphism in C. mauritiana

For the different morphological leaf traits assessed in situ, a weak grouping of trees according to their population of origin was generally observed, except with LA and LDMC, which showed a significant population effect in ANOVA (Analysis of Variance) test (Supplementary Fig. S2). However, the inter-population variability remained difficult to interpret even in LA and LDMC since it was not significantly correlated with the climatic parameters in the population sites (data not shown). To evaluate genetically based variation, inter-population variability was then assessed in the accessions grown in the collection. Substantial variability of morphological traits was found within populations (Supplementary Fig. S3). Multiple and single linear regression tests were performed between both morphological and biochemical traits, and the climatic parameters in the sites where the populations originate, to identify putatively adaptive traits. For the morphological traits, the most significant results in multiple regression (Table�1) were obtained for SLA and LTE (P-value = 4.1 � 10−5 and P-value = 1.2 � 10−6). In particular, variations in temperature (e.g. correlations with Tmin) explained a large part of the variability of both SLA and LTE (54.4% and 52.1%, respectively). Single regression tests (Supplementary Table S2) were also highly significant between SLA and parameters related to temperature (Tmean, Tmax and Tmin) as well as between LTE and Tmin, Tmean or Pmean_JAS. However, no regression tests were yet significant when considering truncated datasets removing all data points from population A, suggesting that the observed correlations were for the most part the result of a single population effect (Supplementary Table S3). Similarly, the significant regression test observed between LDMC and Pmean (P-value= 2.9 � 10−2) was not confirmed when considering a truncated data set without the C population data points. Regarding biochemical traits (Table�1), Alk was the biochemical trait the most significantly associated with variations in climatic factors (P-value = 2.3 � 10−7 and R� = 78.31%). By contrast, variations in wax load were weakly explained by climatic parameters (R� = 51.23% in the global regression model). More specifically, Alk was positively correlated with variations in temperature (e.g. Tmin and Tmean), precipitation (Pmean_JAS) and evaporative demand (PET) (Supplementary Table S2, Supplementary Fig. S4). All these regression tests were still significant when considering truncated data sets that included data from only seven of the eight populations, indicating that the observed correlations did not result from a single population effect (Supplementary Table S3).

Table 1

Multiple linear regression tests between the phenotypic traits assessed in the collection and the climatic parameters at the sites of origin

(a)
Phenotypic trait	Optimized model	Model probability	R�
LA	PET + Pmean	5.7e−03	38.87
SLA	Tmin + Pmean_JAS + PET + Pmean	4.1e−05	72.21
LDMC	Pmean + Tmin	4.3e−04	52.16
LTE	Tmin + Pmean + Pmean_JAS + PET	1.2e−06	77.13
Wax.load	Pmean + PET + Tmax	3.0e−03	51.23
FA	None	–	–
Alk	Tmean + Pmean	2.3e−07	78.31

(a) Global results from optimized regression model (step function). (b) Detailed contribution of each factor included in the ANOVA model.

Significance codes: P-value <0.001: ***; 0.001 < P-value < 0.01: **; 0.01 < P-value < 0.05: *.

Search for cuticle biosynthesis-related genes in C. mauritiana and regression analysis

The biosynthesis of leaf cuticular wax and cutin was then the subject of further investigation. The correlations between the expression profiles of 100 genes involved in cuticule biosynthesis and Alk, the phenotypic trait most closely correlated with climatic parameters, were examined. The expression levels of 10 out of the 100 genes analyzed were significantly correlated with Alk (P-value <0.01). These genes were involved at different steps of wax and cutin biosynthesis (Fig.�4, Supplementary Table S4). In particular, the expression of CER1, which codes for the enzyme catalyzing the final step of alkane biosynthesis, was positively correlated with Alk. By contrast, the expression of CER9, a repressor of alkane biosynthesis, was negatively associated with Alk. Genes coding for LACS and ACT enzymes and ABC transporters were also significantly correlated with Alk, either positively or negatively. Finally, a gene coding for the BODYGUARD protein, which is possibly involved in cutin biosynthesis, was positively correlated with Alk.

Fig. 4

Major steps of cuticular components biosynthesis. The genes whose level of expression was positively associated with the percentage of Alk are in orange boxes. The genes whose level of expression was negatively associated with Alk are in blue boxes. ABC: ,ATP-binding cassette; ACT, acyl-CoA thioesterase; BDG, BODYGUARD; CER, ECERIFERUM; FAH, fatty acyl omega-hydroxylase; FAIH, fatty acyl in-chain hydroxylase; HFADH, omega-hydroxy fatty acyl dehydrogenase; LACS, long-chain acyl-CoA synthetase; OFADH, omega-oxo fatty acyl dehydrogenase; sn-2 GPAT, sn-2 glycerol-3-phosphate acyltransferase.

Analysis of landscape genomics

Since LFMM (Latent Factor Mixed Models) method tests for gene–environment associations while controlling for known population structure (including isolation-by-distance patterns) and demographic history, it was retained for landscape genomics analysis in C. mauritiana. A number of K = 3 latent factors were selected to account for known neutral pattern in this species. Since we were interested in the biological interpretation of the influence of environmental factors on genotypes, we preferentially considered individual climatic parameters showing low correlation among them, instead of using principal components for LFMM analysis. Based on correlations between climatic parameters (Supplementary Table S1, Fig.�1), Trange, Pmean_JAS and Tmin were retained for analysis. A total of 18,321 RNA-seq-derived SNPs obtained across the eight populations sampled in the collection (23 individuals) and filtered for quality criteria were considered in landscape genomics using LFMM analysis.

Running LFMM led to a total of 21 and 17 loci significantly correlated with Pmean_JAS and Tmin, respectively (Supplementary Table S5). Fourteen loci were common to both lists of candidate loci. Expert annotation of the genes underlying significant loci revealed that genes involved in cell wall synthesis and more specifically cell wall permeability, as well as genes involved in abiotic stress response, were particularly well represented. Regarding Trange, LFMM runs did not converged.

Testing candidate SNP loci on an extended panel

To validate the association between three of the most significant SNP loci (located in Cc02_g09660, Cc02_g17520 and Cc02_g07870 genes) and climatic gradients of temperature and precipitation (i.e. Tmin and Pmean_JAS), a wide range of individuals were genotyped by amplicon sequencing. In addition to the accessions used in RNA-seq experiment, 51 individuals across 17 populations covering the spatial distribution of the species (Garot ) were genotyped. The SNPs located on genes Cc02_g07870 and Cc02_g17520 showed differential allelic states according to prevailing climatic conditions, especially regarding Pmean_JAS (Fig.�5). Results of ANOVA and Tukey’s HSD (honestly significant difference) tests supported this result: the mean Pmean_JAS of homozygous genotypes for alternative alleles differed significantly, arguing for an adaptive role of both Cc02_g07870- and Cc02_g17520-related loci at the homozygous state (Table�2, Fig.�5). Allelic frequencies at the heterozygous state were not clearly associated with distinct climatic conditions. Finally, the results of ANOVA were not significant regarding the genotypes at the SNP derived from Cc02_g09660 and thus failed to confirm its role in C. mauritiana adaptation.

Fig. 5

Distribution of the allelic states of three candidate SNPs retained after LFMM analysis, according to the Pmean_JAS and the Tmin. Distribution of allelic states of the candidate SNP supported by the Cc02_g07870 gene according to Pmean_JAS and Tmin is shown in (a) and (b), respectively. Distribution of allelic states of the candidate SNP supported by Cc02_g17520 gene according to Pmean_JAS and Tmin is shown in (c) and (d). Distribution of allelic states of the candidate SNP supported by Cc02_g09660 gene according to Pmean_JAS and Tmin is shown in (e) and (f), respectively. For (a)–(f), the lowercase letters in each figure indicate statistical differences based on Tukey’s HSD test.

Table 2

Results of ANOVA between climatic factors and genotypes of three candidate SNPs retained after LFMM analysis

Climatic factor	CDS	P-value
Pmean_JAS	Cc02_g07870	7.016e−05
Pmean_JAS	Cc02_g17520	0.00471
Pmean_JAS	Cc02_g09660	0.4243
Tmin	Cc02_g07870	0.927
Tmin	Cc02_g17520	0.2531
Tmin	Cc02_g09660	0.9499

Candidate SNPs are designated by the name of the associated CDS. Genotypes were encoded by three modalities: two homozygous states and one heterozygous state. P-values <0.01 are shown in bold.

Discussion

Linking different types of data and analyses makes it possible to contrast results and to refine current hypotheses to provide a better understanding of forest tree adaptation (Stinchcombe and Hoekstra 2008, Sork et al. 2013, Savolainen et al. 2013, De Kort et al. 2014). In the present work, both classical (i.e. common garden experiment) and modern (i.e. transcriptomics and landscape genomics) approaches were successfully combined to uncover processes of local adaptation in a tree endemic to Reunion Island (C. mauritiana). This multi-approach strategy revealed that leaf characteristics appeared to play an important role in the process of local adaptation of the species to contrasted climatic conditions.

Environmentally related differences in leaf morphology are generally considered to reflect an adaptive evolution of plants. For example, plant community-mean leaf size usually decreases with increasing site aridity, arguing for an evolution toward lower respiration and transpiration costs (Givnish 1978, Wright , Ordo�ez ). However, the respective contribution of phenotypic plasticity and local adaptation in controlling variation in leaf morphology often remains unclear. Although substantial variation in leaf morphological traits was observed among populations of C. mauritiana sampled in situ, this diversity was not correlated with climatic factors, thus limiting its interpretation in terms of adaptation. Several factors, including heterogeneity of micro-environment, tree age and pressure of pathogens, can blur potential adaptive traits when assessed in wild populations, thus justifying the need for studies of local adaptation in a controlled environment (de Villemereuil et al. 2016). As the C. mauritiana accessions conserved in the collection have grown in a similar environment since the plantlet stage, the collection can be considered as a common garden experiment. In this context, variation in phenotypic traits among populations is assumed to mainly originate from heritable factors (Gonzalez-Martinez et al. 2006, de Villemereuil et al. 2016). The variation in SLA and LTE among the populations of C. mauritiana sampled in the collection was the most significantly correlated with climatic parameters at their site of origin. Although a trend can be suggested between climate and SLA as well as LTE, the observed environmental gradient for these traits among populations is mainly due to a single population, which may have evolved locally adapted leaf anatomical features to contrasting environmental conditions compared to the other populations. In addition, the low number of populations assessed may also have limited the detection of environmental gradients driving adaptive variation in leaf traits. A study involving a larger number of populations from contrasting habitats would be necessary to confirm or refute these preliminary observations. In addition, a thorough study of leaf anatomy and key physiological processes (i.e. respiration, transpiration and photosynthesis) would be necessary to investigate the ecological significance of variations in leaf traits in C. mauritiana.

By analyzing the composition of cuticular wax among populations conserved in the collection and performing regression tests with climatic parameters, we showed that the relative proportion of alkanes in the wax was the biochemical trait most likely involved in C. mauritiana adaptation to climatic conditions in Reunion Island. Alkane biosynthesis has previously been reported to play a key role in plant response to drought stress. For instance, in Arabidopsis and cotton, plants undergoing water deficit showed an increase in leaf wax load with alkanes accounting for the most part of wax accumulation (Bondada et al. 1996, Kosma et al. 2009). Changes in the proportion of cuticular wax alkane in response to variations in temperature have also been reported in Brassica oleacera (Baker 1974). Changes in the proportion of alkanes in leaf cuticular wax therefore appeared to be a reliable process of local adaptation in C. mauritiana. Leaf wax load appeared to be poorly correlated with climatic variables in C. mauritiana, suggesting that, in this species, adaptation to climatic conditions relies more on a modification of wax composition than on an increase in wax production. Furthermore, the negative correlation observed between the proportions of primary alcohols with that of alkanes suggests that the acyl-reduction and decarbonylation pathways compete in C. mauritiana.

We next used a straightforward approach focusing on the expression of genes involved in the biosynthesis of wax and cutin. The CER1 transcript level was the most significantly correlated with the proportion of alkanes in C. mauritiana leaf wax. In A. thaliana, CER1 is specifically expressed in the epidermis of aerial organs, coexpressed with other genes of the alkane-forming pathway, induced by osmotic stresses and regulated by abscisic acid (Kosma et al. 2009, Bourdenx et al. 2011). CER1 interacts with both CER3 and CYTB5 in the endoplasmic reticulum to catalyze the redox-dependent synthesis of alkanes from very-long-chain acyl-CoAs (Bernard et al. 2012). We also observed that the expression of CER9, which encodes a negative regulator of alkane biosynthesis (L� et al. 2012), was negatively associated with the proportion of alkanes in C. mauritiana.

The levels of expression of genes coding for long-chain acyl-CoA synthetase (LACS), acyl-CoA thioesterase (ACT) and ABC transporters, which are involved in upstream and downstream steps of alkane synthesis, respectively, were also significantly associated with the proportion of alkanes in C. mauritiana cuticular wax. LACS are enzymes that activate FAs released by plastids to CoA-thioesters, prior to wax biosynthesis (Borisjuk et al. 2014). LACS enzymes form a multigenetic family and coexpression of Arabidopsis LACS1 together with the alkane-forming complex CER1–CER3–Cytb5 enhanced alkane production in yeast (Bernard et al. 2012). The positive and negative associations obtained for the expression of LACS7, LACS8 and LACS9 are nevertheless difficult to interpret. It could be that some LACS function in the reverse reaction, hydrolyzing acyl-CoAs thereby producing the free FAs present in wax, together with ACT, which catalyzes the hydrolysis of VLCFA-CoA. ABC transporters are involved in the export of both wax and cutin precursors to the extracellular space (Panikashvili et al. 2007, McFarlane et al. 2010, Bessire et al. 2011). In the same way, ABC transporters form a multigenic family and the specificity of ABC proteins for the translocation of wax precursors remains poorly documented. The positive and negative associations obtained for the expression of ABC1, ABC2 and ABC3 may reflect substrate specificity of the transporters.

The BODYGUARD gene is most probably involved in the extracellular polymerization of cutin precursors in Arabidopsis (Jakobson et al. 2016). In the present study, BODYGUARD expression was positively correlated with the proportion of alkanes in C. mauritiana, finally suggesting that cutin synthesis is simultaneously modulated in association with the modulation of wax biosynthesis.

All these results consistently showed that the biosynthesis of cuticular components, and more specifically the biosynthesis of alkanes, is likely influenced by the climatic conditions at the site of origin of the population concerned. The fact these results were obtained in populations sampled in similar environmental conditions (i.e. in the collection) suggests that the modulation of cuticle composition in C. mauritiana relies on heritable factors. An increase in alkane content has been reported to be involved in plant protection against drought (Bondada et al. 1996, Kim et al. 2007, Kosma et al. 2009, Li et al. 2019). The positive correlation between Alk and both temperature parameters and PET was therefore expected, whereas the positive correlation between Alk and precipitation parameters was not. Reflecting evaporative demand, PET was predictably correlated with Tmin (|r| = 0.80) but less correlated with rainfall parameters (|r| =0.49 with Pmean_JAS and |r|=0.17 with Pmean) implying that well-watered sites can also suffer from atmospheric drought (large vapor pressure deficit and/or high air temperature) in Reunion Island. Thus, the variation in alkane content in C. mauritiana leaves likely reflects local adaptation to warm climates associated with high evaporative demand. The sampling scheme will now be extended to in situ populations, and more specifically to populations encompassing gradients of PET to check the role of alkanes in C. mauritiana adaptation.

Using a landscape genomic approach, we independently searched for additional processes of local adaptation in C. mauritiana. For the purpose of this analysis, SNPs located in coding regions were obtained from RNA-seq experiments conducted in populations sampled in the collection. In addition to being a cost-effective way to identify large numbers of SNPs in non-model organisms, RNA-seq provides SNPs in regions, which are commonly targets of selection and likely underlie phenotypic trait variation (Schork et al. 2013). When untranslated regions (UTRs) are considered, SNPs also cover regions potentially involved in the regulation of gene expression (Barrett et al. 2012).

Several candidate SNP loci revealed by the landscape genomics analysis were located on genes involved in cell wall synthesis and remodeling, including lignin and polysaccharide synthesis. In trees, evidence for natural selection acting on genes affecting cellulose and lignin biosynthesis has frequently been detected by genomic analysis (Dillon et al. 2010, Hadjigol 2012, Camus-Kulandaivelu et al. 2014, Thavamanikumar et al. 2014, Mosca et al. 2016). In addition to being often observed in response to mechanical damage, pathogens or herbivores (Miedes et al. 2014, Agrawal and Weber 2015), the modification in the secondary wall composition with hemicellulose and lignin deposition was reported to be one of the first reactions to abiotic stress (Le Gall et al. 2015).

Genes at significant loci involved in secondary wall remodeling likely coded for: the berberine bridge enzyme (BBE)-like 15, a monolignol oxidoreductase that plays a role in monolignol metabolism and lignin formation (Daniel et al. 2015); the blue-copper-binding (BCB) protein, a lignin synthesis regulator (Ezaki et al. 2005, Ji et al. 2015), which content appears to be associated with cold acclimation in Arabidopsis (Ji et al. 2015); and a glycosylphosphatidylinositol-anchored protein (GAP), whose multigenic family likely has a wide range of functions, including cell signaling, adhesion, cell wall remodeling and pathogen response (Borner et al. 2002, Borner et al. 2003). Finally, significant SNP was also found in genes that were recently evidenced to play a central role in the control of leaf size (KIX8, Li et al. 2018) and the homeostasis of reactive oxygen species and subsequent regulation of stomatal closure in guard cells (i.e. autophagy-associated genes ATG7 and VPS34, Yamauchi et al. 2019, Zhang et al. 2019). Testing the candidate SNP by amplicon sequencing on a larger panel of populations will also check their precise degree of reliability. Among the three candidate loci tested in a larger panel, ANOVA tests revealed that only two SNPs (located on Cc02_g07870 and Cc02_g17520) were significant, suggesting that results of landscape genomics based on a small sampling scheme need further validation using a larger panel. Finally, the interpretation of candidate loci for adaptation was difficult in a non-model species due to the lack of information from gene annotation on multigenic families.

Conclusion

Using complementary approaches, our study provides first insights into processes of local adaptation in a forest tree endemic to Reunion Island. The relative alkane content of leaf cuticular waxes was associated with climatic parameters of the site of origin of the studied populations. Moreover, the differential expression of cuticle biosynthetic genes appeared as a way to constitutively adapt cuticular wax composition to prevailing climatic conditions. Finally, the results of landscape genomic analyses highlighted that genes involved in cell wall composition are also under selection. These findings further our understanding of the processes involved in local adaptation in the indigenous forest trees of Reunion Island and provide clues on the importance of organizing the management of tree populations to cope with climate change.

Materials and Methods

Sampling design

Two separate sampling campaigns were conducted: 8 populations of C. mauritiana were sampled in the international Coffea collection (Saint-Pierre, Reunion Island) (Supplementary Table S6, Fig.�6) and 15 C. mauritiana populations were sampled in forests located across Reunion Island, covering most of the environmental distribution of the species (Supplementary Table S7, Fig.�6). Accessions in the international Coffea collection are produced from plantlets that were collected in situ 10 years ago and transplanted. Thus, the individuals grew in similar climatic conditions until sampling. Leaves of three individual trees per population were collected in collection at the same developmental stage (third leaf from the end of the plagiotropic branch). A total of 12 leaves per tree were used to measure morphological traits and four leaves per tree were used for wax extraction. Populations of C. mauritiana sampled in situ were only used to assess morphological traits. Around seven individuals were sampled per population and eight mature leaves were collected per tree.

Fig. 6

Provenance of C. mauritiana populations sampled in the collection and in situ. The site of origin of populations of accessions sampled in the collection is indicated by red dots. The sites of accessions sampled in situ are indicated by black dots. Populations sampled in both the collection and in situ are indicated by red dots with a black outline. Population names corresponding to population codes are listed in Supplementary Tables S1 and S2. The location of the collection is indicated by a blue triangle.

Measurement of leaf morphological traits

Digital images of leaves collected in both in situ and ex situ populations were produced by scanning. The length (L), maximum width (W) and LA were then determined using ImageJ 1.x software. Leaf fresh mass (FM) was measured directly using a precision balance. After drying in a heat chamber at 50 �C for 48 hours, the leaves were weighed to obtain their dry mass (DM). SLA, LDMC and LTE were calculated according to Vile et al. (2005). Correlations and collinearity between the phenotypic traits (i.e. L, W, LA, FM, DM, SLA, LDMC and LTE) were assessed by computing Pearson’s correlation coefficients and a PCA in R 3.5.2. Non-redundant phenotypic traits were selected among traits that showed correlations with |r| > 0.80 and collinearity in PCA.

Cuticular wax analysis

Eight disks of constant area were die-cut in each leaf sampled in the collection using a 10-mm diameter punch. Cuticular waxes were extracted by immersing the leaf disks for 30 s in chloroform containing 20 �g n-docosane (C22 alkane). Extracts were dried under N2 gas and sent to the CNRS Lipidomics Platform in Bordeaux, France. Wax extracts were then derivatized and analyzed by Gas chromatography-mass spectrometry (GC–MS) to identify the compounds and by GC-FID to quantify the compounds as described in Bourdenx et al. (2011). Leaves of the same tree were considered as biological replicates for wax analysis. The wax.load was measured using n-docosane as internal standard and is expressed in �g per unit of LA.

Climatic data

For each site of origin considered in this study (Fig.�6), a three-decade set of Climate Normals was retrieved from the M�t�o-France R�union database that is fed by a dense network of meteorological stations. The following annual climatic parameters were obtained from climate normals dataset: Tmin, Tmax, Tmean, Trange, Pmean, Pmean_JAS and PET. Pearson’s correlation coefficients were calculated and a PCA was performed in R 3.5.2 to examine possible correlations and to reveal redundant factors (i.e. factors that were correlated at |r| > 0.8 and which showed collinearity in PCA).

Analysis of leaf phenotypic traits

Mean values per tree were used for both leaf morphology and leaf wax composition. Boxplots were performed using the ggplot2 R package (Wickham 2016). Multiple linear regressions were performed between the mean leaf phenotypic traits and the climatic parameters at the site of origin in R 3.5.2 using the step function to optimize the choice of climatic parameters to be included in the model. Simple linear regressions were also performed between population-average foliar phenotypic traits and each climatic parameter individually. To evaluate the robustness of the observed correlations and to rule out the possibility of artifacts due to a single population, each regression model showing P-value <0.01 was reanalyzed, considering truncated datasets removing all data points from one of the eight populations. The eight possible combinations of seven populations were considered. Regarding morphological traits assessed in situ, an ANOVA test was performed in R 3.5.2 to test for the effect of population originating site, followed by Tukey’s HSD post hoc analysis.

RNA sequencing data

Raw data were retrieved from the European Nucleotide Archive no PRJEB31715 (Garot et al. 2019b). Briefly, leaves at the same developmental stage (third leaf from the end of the plagiotropic branch) were sampled in the collection, using the same sampling design as for phenotypic measurements (Supplementary Table S6, Fig.�6), and were immediately frozen in liquid nitrogen. Frozen leaves were then stored at −80�C until RNA extraction. Total RNA extraction, cDNA library construction and single-end sequencing using HiSeq 2500 (Illumina, 100 nt) were performed as described in Garot et al. (2019b). Due to technical difficulties encountered when extracting RNA from a sample, sequencing data were only available for 23 of the 24 accessions. Adapters were trimmed using the cutadapt command (Martin 2011). For the purpose of gene expression analysis, the clean reads were mapped onto the C. canephora reference transcriptome (www.coffee-genome.fr) containing 25,574 coding DNA sequence using BWA-MEM (Li 2013) with the default parameter set, while for the purpose of SNP calling, five prime untranslated (5′-UTR) and three prime untranslated (3′-UTR) regions were considered in addition to CDS in the mapping step. A table of counts of the mapped reads by CDS was obtained using Samtools idxstats (Li et al. 2009). Read counts were normalized using DESeq2 method (Love et al. 2014).

A list of genes involved in the biosynthesis of leaf cuticular wax and cutin in Arabidopsis thaliana was obtained from Li-Beisson et al. (2013). The sequences of these genes were retrieved from the Arabidopsis Information Resource (TAIR) database (www.arabidopsis.org). The X-blast algorithm was used against the C. canephora gene database (www.coffee-genome.fr; Denoeud et al. 2014) to find the best blast mutual hit between Arabidopsis- and C. canephora-predicted peptides. In this way, a total of 176 orthologous genes potentially involved in the biosynthesis of leaf cuticular wax and cutin in coffee were obtained. The expression level of these orthologous genes was retrieved in C. mauritiana. Genes with total read counts <460 (corresponding to a mean count of 20 reads in 23 samples) were filtered out, leaving a total of 100 orthologous genes significantly expressed in C. mauritiana leaves. For each selected gene, linear regressions were performed between the normalized read counts (DESeq normalization) and biochemical traits in R 3.5.2. Results of regression tests were considered as significant at a P-value of <0.01.

SNP filtering

SNP calling was performed according to Garot et al. (2019b) using the Genome Analysis Toolkit pipeline (De Pristo et al. 2011). The following filters were applied to the SNP dataset: SNPs with <10 reads (DP < 10) per individual and Genotype Quality <30 (GQ < 30) per individual were filtered out using BCFtools v.1.3 in the SAMtools package; all SNPs with missing values (not available, NAs), low allelic frequencies (minor allele frequency, MAF < 0.06) and that applied to only one population were removed using VCFtools v. 0.1.13 (Danecek et al. 2011). Using an in-house R script, genotype data were filtered for each gene to only retain one SNP per haplotype, therefore limiting linkage disequilibrium among SNPs. A total of 18,321 SNPs were finally retained.

Testing for associations using LFMM

To detect a possible association between SNP allele frequencies in C. mauritiana accessions and climatic parameters at the site of origin, the latent factor mixed model LFMM v. 2.0 (Caye et al. 2019) was used. Unlike other regression methods, LFMM estimates latent factors (accounting for neutral processes) and regression coefficients simultaneously. This process limits the inclusion of adaptive markers in the ‘neutral set’, which is then used to account for neutral evolutionary processes when estimating (random) environmental effects in regression. Thus, LFMM process reduces the risk of false negatives in test (Frichot et al. 2013). The structure of genetic diversity was assessed using sNMF implemented in LEA package v. 2.0.0 (Frichot and Fran�ois 2015) and the number of latent factors was set to the most reliable number of genetic clusters, as recommended by Frichot et al. (2013), here corresponding to K = 3 (Supplementary Fig. S5). Regularized least squares estimates for latent factor mixed models were calculated using the ridge penalty method with default parameters. Independent runs produced the same results, so only the result of one run is shown. The distribution of P-values was checked to assess the quality of the model (Supplementary Fig. S6). To account for multiple testing, FDR threshold for LFMM outliers was determined using the method of Storey and Tibshirani (2003). The climatic parameters that were previously selected based on PCA results and correlation analysis were used for LFMM analysis.

Annotation

The genes underlying significant SNP loci in LFMM analysis were manually curated. The corresponding protein sequences were compared to known A. thaliana protein sequences by blasting (P-blast) the C. canephora amino acid sequences on TAIR (TAIR10 release). For each sequence, only the best mutual hit was retained. The biological function of the protein was assigned according to the information currently available in the literature.

PCR-based assay for SNP genotyping

A PCR-based assay was used to genotype three candidate SNPs revealed as significant in the LFMM analysis. Genomic DNA for the 23 accessions used in the RNA-seq experiment and for 51 accessions from 17 additional populations (3 accessions per populations) was recovered from the work of Garot et al. (2019a). The latter DNA samples were selected to represent the distribution of genetic diversity of C. mauritiana across Reunion Island. Since C. mauritiana has no sequenced genome, primers for amplifying the region surrounding the three candidate outlier SNP loci were designed on the C. canephora genome using Primer3Plus (Untergasser et al. 2007) (Supplementary Table S8). To guard against polymorphism in targeted primer sequences, the selected primers were also blasted on other available genomes for Coffea species, i.e. on Coffea eugenioides (NCBI, GCA_003713205.1) and Coffea arabica (NCBI, GCA_003713225.1) genomes. SNP detection assays were performed based on Sanger direct sequencing (Eurofins Genomics, Ebersberg, Germany) of DNA amplicons.

Testing the distribution of putatively adaptive alleles

To test for differential distribution of SNP alleles according to climatic gradients, ANOVA was performed between the climatic factors used in the LFMM analysis (values at the site of origin or at the sites where the samples were collected for Pmean_JAS and Tmin) and the genotypes of 74 accessions determined with a PCR-based assay, using the lm function in R 3.5.2 (P-value <0.05). Tukey’s HSD tests were used for the post hoc comparison of means.

Supplementary Data

Supplementary data are available at PCP online.

Funding

The European Union (European Regional Development Fund), R�gion R�union, and the Institut de Recherche pour le D�veloppement (IRD) (CAFEDIV project: ‘Caract�risation et valorisation d’une esp�ce end�mique des for�ts r�unionnaises: le caf�ier marron’; GURDTI/20171191-00002607). E.G. benefited a doctoral grant from the French Ministry of Higher Education, Research and Innovation (University of Montpellier).

Click here for additional data file.