DOI: 10.5061/dryad.6t1g1jx0m
title: Adaptive potential of Coffea canephora from Uganda in response to climate change

description: Understanding vulnerabilities of plant populations to climate change could
 help preserve their biodiversity and reveal new elite parents for future
 breeding programs. To this end, landscape genomics is a useful approach
 for assessing putative adaptations to future climatic conditions,
 especially in long-lived species such as trees. We conducted a population
 genomics study of 207 Coffea canephora trees from seven forests along
 different climate gradients in Uganda. For this, we sequenced 323
 candidate genes involved in key metabolic and defense pathways in coffee.
 Seventy-one SNPs were found to be significantly associated with
 bioclimatic variables, and were thereby considered as putatively adaptive
 loci. These SNPs were linked to key candidate genes, including
 transcription factors, like DREB-like and MYB family genes controlling
 plant responses to abiotic stresses, as well as other genes of
 organoleptic interest, like the DXMT gene involved in caffeine
 biosynthesis and a putative pest repellent. These climate-associated
 genetic markers were used to compute genetic offsets, predicting
 population responses to future climatic conditions based on local climate
 change forecasts. Using these measures of maladaptation to future
 conditions, substantial levels of genetic differentiation between present
 and future diversity were estimated for all populations and scenarios
 considered. The populations from the forests Zoka and Budongo, in the
 northernmost zone of Uganda, appeared to have the lowest genetic offsets
 under all predicted climate change patterns, while populations from
 Kalangala and Mabira, in the Lake Victoria region, exhibited the highest
 genetic offsets. The potential of these findings in terms of ex-situ
 conservation strategies are discussed.
description: Study species and sample selection:<br> Uganda is divided
 into sixteen climate zones based on precipitation patterns as defined by
 Basalirwa (1995), five of which host C. canephora stands. Within these
 five climate zones, 207 georeferenced trees were sampled from seven wild
 forests in 2012 and 2014 by the National Agricultural Research
 Organization (NARO, Uganda) and collaborators of the Institut de Recherche
 pour le Développement (IRD, Montpellier, France). These forests include:
 Budongo (n=65), Itwara (n=23), Kibale (n=19), Kalangala (n=10), Mabira
 (n=25), Malabigambo (n=16) and Zoka (n=49). Populations in Zoka, Budongo,
 Kalangala, Mabira and Malabigambo occurred in distinct climatic envelopes,
 while the climatic envelopes in Itwara tended to overlap those of Kibale
 (Kiwuka et al., 2021). In each targeted forest, leaf samples were
 collected from five sub-sites that were separated by distances of at least
 5 km. Selection of candidate genes and bait
 design:<br> The 323 candidate genes (CGs) selected for the present
 study have been annotated and/or functionally characterized in previous
 studies. They all code for candidate proteins already reported to play
 important roles in central metabolism or in plant responses and adaptation
 to abiotic stress. The CG sequences were retrieved from the whole genome
 assembly of C. canephora (Denoeud et al., 2014) according to the
 annotation available on the Coffee Genome Hub (http://coffee-genome.org/)
 (Dereeper et al., 2015).<br> Probes were designed to cover each CG
 coding region as well as 1 kb upstream and 500 bp downstream flanking
 regions, so as to include putatively regulatory regions. The 120 bp
 MyBaits® probes were designed with 2X tiling and synthesized by MYcroarray
 provider (Ann Arbor, Michigan, USA). A total of 21,306 probes were
 designed. Each candidate probe was BLASTed against the C. canephora genome
 (Denoeud et al., 2014) and filtered based on the manufacturer’s stringent
 criteria (Mariac et al., 2022). Library preparation and
 sequencing:<br> DNA extractions for the 207 samples were performed
 at the IRD facilities from silica-gel dried leaves according to a
 previously described protocol (Mariac et al., 2006). Genomic libraries
 were constructed using the protocols outlined in Rohland & Reich
 (2012) and Mariac et al. (2014). The 207 individual libraries were then
 capture-enriched by pools of 48 libraries using the synthetic RNA MyBaits®
 probes and according to the MYcroarray protocol (Mariac et al., 2022). The
 enriched pools were quantified using real-time PCR and combined in
 equimolar ratios prior to sequencing on one lane of 150 bp paired end
 reads on an Illumina HiSeq 3000 sequencer (GeT-PlaGe Platform, GenoToul,
 Toulouse, France). SNP genotyping, calling and
 filtering:<br> Sequence analysis was performed using scripts
 published by Mariac et al. (2014) and Scarcelli et al. (2016) and also
 available on GitHub (https://github.com/Maillol/demultadapt;
 https://github.com/SouthGreenPlatform/arcad-hts/blob/master/scripts/arcad_hts_2_Filter_Fastq_On_Mean_Quality.pl).<br> The mapping step was carried out using BWA MEM 0.7.5a-r405 (Li & Durbin, 2009) with the default option (-B 4) and the C. canephora assembly (http://coffee-genome.org/coffeacanephora) as reference. SNP calling was done using UnifiedGenotyper in the Genome Analysis Toolkit (GATK v3.6). SNPs located on the selected CG sequences were considered as ‘in-target’ and the other ones as ‘off-target’. Two successive sets of filters were applied to raw SNPs. We first discarded low quality variants according to the quality criteria recommended by GATK, and selected only biallelic SNPs using VCFtools v0.1.13 (Danecek et al., 2011).<br> We applied additional filters for population genetic analyses and for association analyses, i.e. keeping SNPs with no excess of heterozygous genotypes (&lt; 0.8), a minor allele frequency (MAF) greater than 5% and under linkage equilibrium. For the latter filter, SNPs were processed with PLINK 1.90b4 (Purcell et al., 2007) to prune only SNPs in approximate linkage equilibrium based on the pairwise correlation between the SNP genotype counts for 100 bp sliding windows with 10 bp steps (option -indep-pairwise). The SNPs were considered correlated when r2 &gt; 0.5. These filters led to a total of 5,860 SNPs: 4,753 in-target and 1,107 off-target loci. Bioclimatic data:<br> Environmental factors (bioclimatic variables BIO1-19, Table S1) were downloaded from the WorldClim database (http://www.worldclim.org, Fick & Hijmans, 2017) at 30 arc-second resolution (~1 km) for ‘Current conditions ~1960-2000’