Received: 5 December 2018 | Revised: 18 June 2019 | Accepted: 24 June 2019 DOI: 10.1111/eva.12838 ORIGINAL ARTICLE A multiscale approach to detect selection in nonmodel tree species: Widespread adaptation despite population decline in Taxus baccata L Maria Mayol1 | Miquel Riba1,2 | Stephen Cavers3 | Delphine Grivet4,5 | 6 7 8 Lucie Vincenot | Federica Cattonaro | Giovanni G. Vendramin | 1,9 Santiago C. González‐Martínez 1 CREAF, Cerdanyola del Vallès, Spain 2 Abstract Univ. Autònoma Barcelona, Cerdanyola del Vallès, Spain Detecting the molecular basis of local adaptation and identifying selective drivers 3 Centre for Ecology and Hydrology is still challenging in nonmodel species. The use of purely population genetic ap‐ Edinburgh, Penicuik, UK 4 proaches is limited by some characteristics of genetic systems, such as pleiotropy Department of Forest Ecology and Genetics, Forest Research Centre, INIA‐ and polygenic control, and parallel evidence from phenotypic‐based experimental CIFOR, Madrid, Spain comparisons is required. In long‐lived organisms, the detection of selective pressures 5 Sustainable Forest Management Research might also be precluded by evolutionary lag times in response to the environment. Institute, INIA‐University of Valladolid, Madrid, Spain Here, we used the English yew to showcase an example of a multiscale integrative 6 UNIROUEN, IRSTEA, ECODIV, Normandie approach in a nonmodel species with limited plant and genomic resources. We com‐ Univ., Rouen, France 7 bined information from two independent sources, phenotypes in a common envi‐ Instituto di Genomica Applicata (IGA), Udine, Italy ronment and genomic data in natural populations, to investigate the signature of 8 Institute of Biosciences and Bioresources, selection. Growth differences among populations in a common environment, and Division of Florence, National Research Council, Sesto Fiorentino, Italy phenological patterns of both shoot elongation and male strobili maturation, were 9 UMR BIOGECO, INRA, Univ. of Bordeaux, associated with climate clines, providing evidence for local adaptation and guiding us Cestas, France in the selection of populations for genomic analyses. We used information on over Correspondence 25,000 SNPs from c. 1,200 genes to infer the demographic history and to test for Maria Mayol, CREAF, Cerdanyola del Vallès molecular signatures of selection at different levels: SNP, gene, and biological path‐ 08193, Spain. Email: maria.mayol@uab.es way. Our results confirmed an overall demographic history of population decline, but we also found evidence for putative local adaptation at the molecular level. We Funding information Spanish Ministry of Economy and identified or confirmed several candidate genes for positive and negative selection Competitiveness, Grant/Award Number: in forest trees, including the pseudo‐response regulator 7 (PRR7), an essential compo‐ ADAPCON (CGL2011-30182-C02-01/02) and TAXUS (CGL2007-63107⁄ BOS); nent of the circadian clock in plants. In addition, we successfully tested an approach European Union’s Horizon 2020 research to detect polygenic adaptation in biological pathways, allowing us to identify the fla‐ and innovation programme, Grant/Award Number: 676876 (GenTree) vonoid biosynthesis pathway as a candidate stress‐response pathway that deserves further attention in other plants. Finally, our study contributes to the emerging view that explaining contemporary standing genetic variation requires considering adap‐ tation to past climates, especially for long‐lived trees. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2019 The Authors. Evolutionary Applications published by John Wiley & Sons Ltd Evolutionary Applications. 2020;13:143–160.  wileyonlinelibrary.com/journal/eva | 143 144 | MAYOL et al. KEYWORDS adaptation, demographic decline, English yew (Taxus baccata), environmental association, polygenic adaptation, single nucleotide polymorphism 1 | I NTRO D U C TI O N al., 2014), and even fewer in the case of nonmodel tree species with limited available genomic resources. Understanding adaptive genetic responses to climate is crucial for The characteristics of genetic systems, such as pleiotropy and preserving biological diversity and predicting future evolution‐ polygenic control, place limits on the use of purely population ge‐ ary responses, particularly for many forest trees that are keystone netic approaches to infer local adaptation, known as “reverse ecol‐ species within the ecosystem. In the past few years, the increase ogy” (Li, Costello, Holloway, & Hahn, 2008). To successfully use in forest tree genomic resources, coupled with the development of genomic analyses to detect the molecular basis of local adaptation novel analytical approaches, has contributed to understand how and identify selective drivers require support from organism‐based adaptive and neutral processes have shaped genetic variation of tree approaches (Tiffin & Ross‐Ibarra, 2014). Unfortunately, in many species (Holliday et al., 2017). Nevertheless, there are some major cases, genomic analyses targeting adaptation often do not con‐ challenges ahead. On the one hand, from a methodological point of sider relevant additional information on phenotypic traits (Barrett & view, the power to detect loci under selection with commonly used Hoekstra, 2011; Villemereuil, Gaggiotti, Mouterde, & Till‐Bottraud, approaches (e.g., FST‐outlier detection or environmental association) 2016). In forest trees, common garden approaches have been widely is variable depending on the sampling effort and the number of en‐ used to study adaptive evolution in traits such as growth, wood vironmental variables included (Ahrens et al., 2018). On the other properties, and cold‐hardiness, but most experiments have been re‐ hand, and despite the considerable effort made to discern between stricted to species and traits of economic value. Extending the use selective and purely demographic processes, these methods are of experiments to species of noneconomic interest, especially those still prone to high false‐positive rates (Lotterhos & Whitlock, 2014; with contrasting life‐history traits, would improve our understand‐ Villemereuil, Frichot, Bazin, Francois, & Gaggiotti, 2014). Moreover, ing of the multiple factors that influence adaptation in forest trees. accounting for the underlying population structure to reduce false However, performing well‐designed, replicated common garden ex‐ positives may end up removing true adaptive variants when popula‐ periments on nonmodel species, or those not considered economi‐ tion structure and environmental drivers of adaptation are aligned, cally important, is certainly a challenge in view of the general lack particularly when the trait is polygenic and/or the selection signa‐ of funding. This is even more problematic when, due to the species’ ture is weak (Ahrens et al., 2018). Some alternative strategies have distribution and population characteristics, access to plant material been applied to circumvent these problems. For example, in a study is limited. investigating convergent adaptation to climate in distantly related Finally, association of climatic variables with potentially adaptive conifers (Pinus contorta and Picea glauca complex), the authors did molecular variation has usually been restricted to using contempo‐ not correct by population structure but assumed that significant rary environmental data. However, in tree populations, time lags environmental associations shared across species could not have may develop in genotype–environment associations and persist for arisen from purely neutral processes (Yeaman et al., 2016). Similar many generations as a consequence of long generation times and approaches can then be applied to complement conventional meth‐ late reproductive age. Therefore, genetic variation in long‐lived ods of studying adaptation with genomic data. species might not be in equilibrium with current environments but From a biological point of view, an additional and particularly rather be associated with past ones (e.g., Gugger, Ikegami, & Sork, important limitation is that most methods are designed to detect 2013; Epps & Keyghobadi, 2015; Mayol et al., 2015). Exploring the the signature of selection at individual loci. However, many adap‐ association of genetic variability with environment through time is tive traits are probably governed by multiple loci, each contributing essential to dissect the relative contribution of past and present cli‐ only modestly to the phenotype, and adaptation is expected to take mates to the patterns observed. This can be addressed by applying place by selection on standing variation at many loci simultaneously a temporal dimension to population genomic approaches: Loci con‐ (Pritchard, Pickrell, & Coop, 2010). Consequently, in model species ferring adaptation to past environments that are selectively neutral like humans, an increasing number of studies have shifted from fo‐ in modern‐day conditions would be expected to show significant as‐ cusing on single genes to polygenic approaches based on biologi‐ sociations with paleoenvironments, while the inverse would be ex‐ cally meaningful groups of genes that underpin processes such as pected for loci conferring adaptation to current climate (Lafontaine, metabolic pathways (e.g., Daub et al., 2013; Daub, Moretti, Davydov, Napier, Petit, & Hu, 2018). Excoffier, & Robinson‐Rechavi, 2017; Berg & Coop, 2014; Foll, In this study, we used the English yew (Taxus baccata L.) to show‐ Gaggiotti, Daub, Vatsiou, & Excoffier, 2014). To date, however, few case an example of an integrative approach in a nonmodel species studies have applied multilocus approaches to investigate the mech‐ with limited plant and genomic resources, taking into consideration anisms underlying local adaptation in forest trees (but see Csilléry et the issues raised above concerning the study of local adaptation. We MAYOL et al. | 145 first addressed the phenotypic variability and potential adaptive re‐ Koutsodendris et al., 2010) and molecular data, at least for popula‐ sponses of trees from a reduced but climatically heterogeneous part tions from the Iberian Peninsula (Burgarella et al., 2012). of the species range (Iberian Peninsula) growing in a common envi‐ The small populations typical of English yew are predicted ronment. From the results obtained, we designed an independent to have restricted capacity to adapt to environmental change be‐ geographic sampling to enable a rangewide genomic analysis. Taking cause enhanced genetic drift opposes the effect of natural selec‐ advantage of the recently available T. baccata transcriptome (Olsson tion (Wright, 1931) and increases genetic load (i.e., fixed deleterious et al., 2018), we designed a gene capture approach for ~1,200 genes, mutations; Whitlock, 2000). However, some empirical studies have including known candidate genes for climate adaptation in forest suggested that divergent selection on phenotypic traits might also trees. Then, we inferred the demographic history of populations and be enhanced in small, fragmented populations (Willi, Van Buskirk, applied several methods to investigate the signature of selection at Schmid, & Fischer, 2007). Consistent with these results, the signa‐ single‐locus (SNP) and multilocus (gene, biological pathway) levels. ture of ongoing or recently completed selection was found in some Since neutral genetic diversity at the whole distribution level bears of the genes involved in taxol biosynthesis in Iberian yew popula‐ the signature of potential adaptation to present and past environ‐ tions, despite a clear signature of recent (2,000–3,000 generations) ments (Mayol et al., 2015), we used climatic information from the bottlenecks (Burgarella et al., 2012). Our goals were thus (a) to in‐ present, the last glacial maximum and the last interglacial periods to vestigate whether other potentially adaptive traits, such as those evaluate the relative importance of current and past climatic condi‐ related to growth and phenology, show evidence of local adaptation tions on patterns of genetic variation. in Iberian populations of the species, (b) to assess the role of local ad‐ English yew is a slow‐growing and long‐lived dioecious gymno‐ aptation at the rangewide scale in a species characterized by demo‐ sperm presently found in a variety of habitats throughout most of graphic bottlenecks, and (c) to determine the contribution of present the European continent. Despite its wide distribution, it often forms and past environments to the observed patterns of variation and the small and isolated populations, mainly in the Mediterranean area most likely environmental drivers underlying them. (Thomas & Polwart, 2003). Most natural or seminatural populations are not easily accessible, and obtaining large amounts of seed per 2 | M ATE R I A L S A N D M E TH O DS population and female tree is a limiting issue for adequate replication in provenance trials. Furthermore, conditions for seed germination 2.1 | Phenotypic variability in a common garden and storage are highly restrictive, and dormancy is highly variable (Thomas & Polwart, 2003). In contrast to other gymnosperms, gene The analysis of phenotypic variation was carried out in the Valsaín flow among English yew populations is limited, and genetic diversity clonal bank, located in central Spain (Segovia, 40º54′38″N, is highly structured both at the local (Dubreuil et al., 2010) and the 04º00′45″W, 1,136 m a.s.l; Figure 1). This unique clonal T. bac‐ regional scale (Chybicki, Oleksa, & Kowalkowska, 2012; González‐ cata collection has been progressively built since 1992 for ex situ Martínez et al., 2010). Its demographic history seems to be mostly conservation of genetic resources in the Iberian Peninsula. Plant characterized by a strong and continuous reduction of the effec‐ material (genets) was collected in 26 native populations from one‐ tive population size since the last interglacial, as suggested by pa‐ year‐old shoots located on terminal branches. Due to the patchy laeoecological records (e.g., Turner, 2000; de Beaulieu et al., 2001; distribution of the species, populations are easily defined as a F I G U R E 1 Geographic distribution of the 26 natural populations of Taxus baccata planted in the Valsaín clonal bank (see details in Table S1). Colors for populations indicate different bio‐ geographic regions: Atlantic Region (blue dots), Cantabrian Mountains (dark pink dots), Catalan Coastal Range (orange dots), Central System (green dots), Pyrenees (yellow dots), and Sub‐Baetic System (red dots). The location of Valsaín experimental site is indicated by a black star 146 | MAYOL et al. group of reproducing individuals growing along ravines or hillsides (early and mid/late March) on tagged apical shoots and was estimated and readily identified in the landscape. Clones were propagated as the proportion of strobili shedding pollen using pooled data from through stem cuttings in growth chambers and later transferred both samplings. In both cases, we used mixed models to determine outdoors. Plants were planted 2–3 m apart in a row and column the effects of Population, Sex and Year, using genet and individual layout without following any particular randomized design, re‐ within genet as random effects. In the analysis of shoot elongation, placed when dead, and watered if needed during summer. After we also tested for Sex by Year and Population by Year interactions. accounting for ramets failing to establish, the available plant ma‐ Models were conducted on the original scale assuming a Gaussian terial included a total of 254 uneven aged (10–15 year old) clones distribution of residuals. Strobili maturation was assessed in a subsa‐ (individual plants) belonging to 108 genets (1–9 genets per popu‐ mple of 21 populations. Since several individuals and genets failed to lation) distributed across six contrasted geographic regions in the reproduce every year, we performed 2 analyses. In the first, we used Iberian Peninsula (Figure 1, Table S1). 14 populations allowing us to test for the effect of Population, Year, First, we measured shoot growth, which is often considered as a and their interaction. In the second, we tested only for Population fitness proxy in forest trees (e.g., Alía, Chambel, Notivol, Climent, & and Year main effects using all populations. Models of the proportion González‐Martínez, 2014). Since young yew trees tend to grow as a of mature strobili were run assuming binomial random errors. shrub, individual growth was estimated as the volume (diameter2 * All mixed models were performed in R v. 3.4.4 (R Core Team, 2018) length, in mm3) of the leading current‐year shoot located in the two using the lme4 package (Bates, Maechler, Bolker, & Walker, 2015). For dominant (longest) stems in four consecutive years (2009–2012). models with Gaussian errors, we report F tests based on Kenward– We used the average value of both leading shoots as the response Roger's approximation of the degrees of freedom. Reported P‐values variable. According to the number of replicates available, the vari‐ for mixed models with binomial type errors are based on likelihood ability in growth was assessed using a Gaussian mixed model with ratio tests. Comparisons of mean values across groups were performed Population, Sex, Year, and Population by Year interaction as cate‐ using the emmeans package (Length, 2018). Estimated marginal means gorical fixed factors. We used stem length at the beginning of our for each Population were used to determine associations with climate follow‐up study (2009) as a covariate. Genet and individual within variables. Monthly records of temperature (mean, maximum and min‐ genet were treated as random factors. Stem length and shoot vol‐ imum) and precipitation were obtained from the Digital Climatic Atlas ume were log‐transformed after residual analyses. of the Iberian Peninsula (http://openg​is.uab.es/wms/iberi​a/en_index. Second, we investigated the phenology of shoot elongation and htm; Ninyerola, Pons, & Roure, 2005). From these data, we calculated male strobili maturation from 2010 to 2012. Growth phenology was monthly temperature range and averaged seasonal values for climatic assessed by comparing the ratio of shoot length produced before variables: winter (December‐February), spring (March‐May), summer the onset of summer drought (May‐June) to the total shoot length (June‐August), and fall (September‐November). For the same variables attained at the end of the growing season (October‐November). and periods, we also included estimates of calibrated climatic distance Measures of length were taken on the leading current‐year shoot of (Gower's distance) between the location of each source population two dominant stems per individual and averaged before analyses. and the average values at the common garden site recorded during Male strobili maturation was assessed in two consecutive samplings the study period (see Rutter & Fenster, 2007). F I G U R E 2 Location of the European populations of Taxus baccata sampled for the gene capture experiment. Pie charts show the proportion of gene pool membership (K = 2) by population inferred from 25,726 SNPs with faststructure. See Table S2 for details on populations MAYOL et al. | 147 sequencing of captured regions on an Illumina platform. Probes were 2.2 | Sampling for molecular analyses and designed by Roche Tech‐Support in Madison (WI, USA) starting from DNA extraction the selected 1,134 gene sequences (about 1.6 Mbp of target se‐ Phenotypic variability in the common garden suggested a potential quence in total). Library preparation and targeted sequencing were cline of adaptation from cold, inland, continental climates to milder, outsourced to IGA Technology Services (Udine, Italy). Libraries were temperate, and coastal ones (see Results). Based on these results, constructed by using KAPA DNA Library Preparation Kit (Roche, we sampled 12 populations throughout the European distribution NimbleGen Inc.) following the manufacturer's protocol and enrich‐ of English yew spanning a range of variation from cold continental ment performed using NimbleGen solution‐based SeqCap EZ probe to milder temperate conditions (Figure 2, Table S2). Since we aimed libraries kit. Cluster generation, template hybridization, isothermal to test for the effects of both present and past climatic conditions, amplification, linearization, blocking and denaturization, and hybridi‐ populations were sampled from locations where suitable environ‐ zation of the sequencing primers were then performed on Illumina ments for yew persistence were present during the last glacial cBot and flow cell HiSeq SBS V4 50 cycle kit, loaded on HiSeq 2500 maximum (LGM, c. 21,000 year BP) and the last interglacial (LIG, c. Illumina sequencer producing 50 bp single reads. The CASAVA ver‐ 120,000–140,000 year BP), in agreement with Mayol et al. (2015). sion of the Illumina Pipeline 1.8.2 was used for base calling and de‐ Ten individuals were sampled in each population, resulting in a total multiplexing. Adapters were masked using cutadapt (Martin, 2011). sample size of 120 individuals. Total DNA was isolated from 50 to Masked and low‐quality bases were filtered using erne‐filter (Del 100 mg of dry leaf material using the Invisorb® DNA Plant HTS 96 Fabbro, Scalabrin, Morgante, & Giorgi, 2013). Kit (STRATEC Molecular GmbH, Berlin, Germany). The GATK pipeline was used for SNP calling, following GATK best practices for SNP filtering (–filterExpression “MQ0> = 4 && ((MQ0/ (1.0 * DP)) > 0.1) || DP < 10.0 || Q < 50.0 || QD < 1.5 || FS > 60.0”), 2.3 | Selection of genes for sequencing which allowed the initial identification of 52,904 SNPs. SNP data We used two main strategies to select candidate genes. First, we re‐ were further filtered using VCFtools v0.1.14 (Danecek et al., 2011) viewed the literature (up to 2012) focused on the characterization of and vcffilter (part of vcflib C++ library) for average coverage across loci putatively involved in the adaptation of tree species to environ‐ all samples (between 20 and 250) and allele balance (between 0.3 ment (e.g., cold, drought, pathogens, phenology) and obtained a list and 0.7 for the heterozygous calls). Only polymorphic biallelic SNPs of 344 genes of interest (Table S3). Then, we retrieved the reference with fewer than 15% missing data were kept. This stringent filtering sequences from these publications and aligned them to the T. bac‐ reduced the SNP dataset to 36,028 SNPs. This dataset was addi‐ cata transcriptome (Olsson et al., 2018) using BLAST (Altschul, Gish, tionally filtered by removing SNPs that showed a significant hetero‐ Miller, Myers, & Lipman, 1990). Positive results were obtained for zygote excess in ≥4 populations to minimize the impact of putative 273 genes > 450 bp, which were selected for subsequent analyses. paralogous loci on data analysis, leaving a final dataset of 25,726 Second, the transcriptomes of T. baccata (Olsson et al., 2018), SNPs spanning 1,128 genes (six out of the initial 1,134 genes were T. cuspidata (Lee et al., 2010), T. wallichiana (Hao, Ge, Xiao, Zhang, & discarded because of low sequence quality). Yang, 2011), and Picea sitchensis (European Bioinformatics Institute, TSA project GACG01000000 data at http://www.ebi.ac.uk/) were 2.5 | Genetic diversity, population used to compute Ka/Ks and detect additional potential candidates structure, and demography under moderate/high positive evolution (Ka/Ks > 0.5 for comparisons with T. cuspidata and T. wallichiana; and Ka/Ks > 0.3 for those with Overall nucleotide diversity (π), the efficacy of selection evaluated P. sitchensis). From this approach, we obtained 345 candidate genes as the ratio of nonsynonymous to synonymous nucleotide diversity with > 450 bp length and functional annotation available, so the final (πa/πs), and neutrality tests (Tajima's D, and Fu & Li's D* and F*) were number of candidate genes obtained from literature and sequence computed using mstatspop (https​://bioin​forma​tics.crage​nomica.es/ divergence ratios was 618 (797,452 bp). numge​nomic​s/peopl​e/sebas/​softw​are/softw​are.html) on concat‐ Finally, we completed our design including 516 genes with enated sequence files, both considering all populations together > 450 bp length for which functional annotation was available and separately for distinct geographic regions (Bosnia‐Herzegovina, (753,183 bp), despite not being selected as candidates initially by Iberian Peninsula, Slovakia, Switzerland and United Kingdom, Table any of the two approaches. Thus, the final number of genes included S2). Average inbreeding coefficients (F) and pairwise relatedness in the study was 1,134 (1,550,635 bp). (unadjusted Ajk statistic, Yang et al., 2010) were computed using VCFtools v0.1.14 (Danecek et al., 2011). Additionally, relatedness between pairs of individuals within each population was computed 2.4 | Probe design, library construction, and to discard the potential inclusion of replicated individuals (clones) in targeted sequencing the samples. Selected genes were sequenced in the 120 yew trees by a se‐ Population genetic structure was analyzed using the program quence capture approach using a custom SeqCap EZ design (Roche faststructure (Raj, Stephens, & Pritchard, 2014), from K = 1 (i.e., NimbleGen Inc, Basel, Switzerland) followed by next‐generation no structure) to K = 12 (the number of localities in our sample). 148 | MAYOL et al. Runs were repeated three times, and averaged Q‐values (i.e., the frequency of SNPs related to climate per gene to be P = ∑iai/ni individual assignment probability to each of the K groups) were (summation over i genes), and calculating P separately for each used to draw pie charts in Figure 2. In a previous study, Burgarella climatic variable. Following Yeaman et al. (2016), genes with zero et al. (2012) showed strong demographic decline (based on nuSSRs) outliers were excluded from the calculation of P. Finally, to control in yew populations from the Iberian Peninsula. To confirm or for false positives and produce a more robust test, we combined refute these results at a wider geographic range, we assessed effec‐ the results of this approach with those obtained with the SNP tive population size over time using the stairway plot approach on analyses and retained only those genes that also had at least two unphased SNP data (version 2 beta; Liu & Fu, 2015). We assumed different SNPs that were significantly associated with the same a per‐generation mutation rate of 1.4 × 10–8 considering per‐year climatic variable in bayescenv and/or baypass . Candidate genes were mutation rates for conifers in Chen, Uebbing, et al. (2012) and a annotated using BLASTX searches against the Picea abies genome generation time for English yew of 50 years (Thomas & Polwart, at ConGenIE public database (http://www.conge​nie.org; Nystedt 2003). Median estimates of the population size and confidence et al., 2013), using a minimum threshold of 50% of sequence iden‐ intervals were estimated with the built‐in bootstrap function using tity and E‐value < 1E‐50. The functional annotation of candidate 200 replicates of the input file. genes was then based on the associated A. thaliana description for each of the best hits provided with the Norway spruce reference genome sequence. 2.6 | Detecting signatures of selection at SNP, Third, we used a complementary, independent, procedure to gene, and pathway level investigate the potential joint effect of selection acting on multiple We combined several approaches to investigate the contribution loci (polygenic adaptation) instead of on single genes. This procedure of natural selection to genetic adaptation in English yew. First, was, in turn, independent of the environmental variables chosen, al‐ we searched for signals of selection at the single‐locus (SNP) level lowing the detection of potential signatures of selection in response using two environmental association approaches that take into ac‐ to other drivers apart from climate. For this purpose, we used KEGG count population genetic structure to avoid false positives: bayes- Automatic Annotation Server to map our 1,128 genes and retained cenv (Villemereuil & Gaggiotti, 2015) and baypass (Gautier, 2015). 75 pathways that included at least five genes. We refined this path‐ For each population, we downloaded the climatic information avail‐ way list by retaining only known pathways in plants, most of them able at the WorldClim database (Hijmans, Cameron, Parra, Jones, & having an Arabidopsis thaliana reference map. The final dataset thus Jarvis, 2005) for the present (PRE, c. 1950–2000), the last glacial consisted of 390 distinct genes that were mapped on 60 pathways, maximum (LGM, c. 21,000 year BP) and the last interglacial (LIG, c. of which eight were plant‐specific. Then, for each sequence‐based 120,000–140,000 year BP) periods. We obtained data for eight cli‐ statistic obtained with mstatspop (π, πa/πs, Tajima's D, and Fu & Li's D* matic variables for each period: annual mean temperature, tempera‐ and F*), we performed a gene set enrichment approach as in Daub et ture seasonality (standard deviation), maximum temperature of the al. (2013), using polysel gene set enrichment pipeline freely available warmest month, minimum temperature of the coldest month, annual from https​://github.com/CMPG/polysel. Briefly, for each statistic, we precipitation, precipitation of the wettest month, precipitation of calculated the SUMSTAT score (Tintle, Borchers, Brown, & Bekmetjev, the driest month, and precipitation seasonality (coefficient of vari‐ 2009), which is simply the sum of the values for that statistic of all ation). bayescenv was run for each of the standardized climatic vari‐ genes in a given pathway, and compared it to a null normal distribution ables with general default parameters. After 20 pilot runs of 2,000 of 100,000 random sets with the same size as the original pathway. iterations and a burn‐in of 50,000 iterations, 5,000 MCMC samples The null distribution was created with a sequential random sam‐ were taken with 10 steps between each sample. For baypass runs, pling method (Ahrens & Dieter, 1985). Following Daub et al. (2013), standardized variables were run using default parameters under the a “pruning” approach was applied to avoid any redundancy between standard model (STD in Gautier, 2015). We retained all significant gene sets (see also polysel manual). Briefly, the genes belonging to the SNP–climate associations with Q‐value < 0.05 (bayescenv) and/or highest scoring gene set were removed, iteratively, from the remain‐ Bayes Factor > 30 (baypass). ing gene sets and the testing procedure was re‐run until no gene sets Second, we searched for genes with a high proportion of SNPs large enough to be tested were left. When following this procedure, being correlated to climatic variables using a similar approach to tests are not independent and biased toward low P‐values (i.e., only the method proposed by Yeaman et al. (2016). For each SNP, we the higher scoring gene sets remain after pruning), and thus, we esti‐ calculated Spearman's rank correlations between allele frequency mated the false discovery rate (FDR) empirically using 300 iterations. and the 24 variables (8 climatic variables × 3 periods), and retained all SNPs with significant correlations at p < 0.01. For each gene, 3 | R E S U LT S we then counted the total number of SNPs (n) and the number of SNPs with significant climate correlations (a). Then, to identify 3.1 | Adaptive variation at the phenotypic level candidate genes potentially under selection, we compared the number of SNPs related to climate in each gene to the 0.9999 The analysis of growth showed significant differences between quantile of its binomial expectation, considering the expected sexes (F1, 66.03 = 4.65; p < 0.05). Males showed higher annual growth MAYOL et al. | 149 than females (Figure S1). We also found differences among popula‐ tions across years (Population by Year interaction: F75, 650.02 = 1.85; p < 0.001). However, relative differences among populations were quite consistent across years (Figure S2). We also found signifi‐ cant negative correlations between mean population growth and calibrated mean and maximum winter temperatures (Pearson's r = −0.58 and −0.49, respectively; p < 0.05; Figure 3a), suggesting local adaptation for growth to temperatures experienced during the colder seasons. Removing the populations with the lowest level of replication (those with less than three genets) produced the same re‐ sults and improved the relationship between growth and calibrated environmental variables (e.g., r = −0.74 and p < 0.001 for mean win‐ ter temperature). Significant differences among populations were also found in phenological patterns of shoot elongation and maturation of male strobili. The proportion of late shoot growth (from early summer to late fall) varied among populations and between sexes, but such dif‐ ferences were also dependent on the year considered (Population by Year interaction: F50, 389.97 = 1.89; p < 0.001; Sex by Year interaction: F2, 391.32 = 4.49; p < 0.05). However, differences among populations and sexes were highly consistent across years (Figures S3 and S4, respectively). Compared to males, higher values of late shoot growth in females suggest more protracted and slower growth in females. Patterns of variation in shoot growth phenology were significantly associated with temperature range during summer and fall (Pearson's r = −0.56 and −0.51, respectively; p < 0.05; Figure 3b). This relation‐ ship suggests that yew trees from more continental sites or located at higher elevations, with largest diurnal temperature range, started growth earlier and/or grew faster. The analysis of the variability in male strobili maturation (using a subset of populations, see Materials and Methods) showed a significant interaction between population and year (χ 2 = 476.64, df = 26, p < 0.001), although the overall pattern of differentiation among populations was also similar across years (Figure S5). Differences among years and populations in strobili phe‐ nology were also significant when considering all available popula‐ tions (Year: χ2 = 1,770.5, df = 2, p < 0.001; Population: χ2 = 44.64, df = 20, p < 0.01). We also found a significant negative correlation between strobili maturation and temperature, particularly mean and maximum temperature during winter and spring (Pearson's r= −0.64 to −0.61, p < 0.01 in all cases; Figure 3c), with male strobili maturing earlier in those genets from colder populations. F I G U R E 3 Phenotypic variation measured in the common garden (Valsaín clonal bank) and its correlation with climate. (a) Correlation 3.2 | Genetic diversity, population between mean population shoot growth (volume, mm3) and calibrated structure, and demography winter temperature (Gower's distance); (b) correlation between the mean proportion of late (summer‐fall) shoot growth (in length) and Overall synonymous nucleotide diversity (πs = 0.00897) and efficacy mean temperature range in summer; (c) correlation between male of selection (πa/πs = 0.308; Table 1) were similar to those found in strobili maturation (estimated as the proportion of open strobili other outcrossing long‐term perennial plants (reviewed in Chen, shedding pollen relative to the total number of strobili) during the reproduction peak and spring maximum temperature. Models are based Glemin, & Lascoux, 2017). At the regional level, Slovakia and Bosnia‐ on estimated marginal means, see Materials and Methods for details Herzegovina (Eastern gene pool, see below and Mayol et al., 2015) showed a lower proportion of polymorphic sites, probably because that is, Bosnia‐Herzegovina and the Iberian Peninsula. None of the the smaller number of individuals sampled in these regions (Table 1). individuals sampled was suspicious to be replicated due to clonal‐ Nucleotide diversity was consistently lower at southern latitudes, ity, since high values of pairwise relatedness (>0.6) within each 150 | MAYOL et al. population were restricted to the comparisons of an individual with TA B L E 1 Genetic diversity, neutrality test, and inbreeding statistics. Main statistics are given at the species level and separately for the different Taxus baccata geographic regions included Related −0.008 −0.088 −0.030 −0.030 −0.087 −0.023 itself (Table S4). in this study. πa/πs: efficacy of selection as evaluated by the ratio of nonsynonymous to synonymous π (Tajima, 1989); TajD: Tajima's D (Tajima, 1989); FuD*: Fu & Li's D* (Fu & Li, 1993); FuF*: Fu & Li's F* (Fu & Li, 1993); F: average inbreeding coefficient; Related: average pairwise relatedness statistic (unadjusted Ajk) based on the method of Yang et al. (2010); values around zero Population genetic structure showed an optimal partition into two genetic pools located at the eastern (Bosnia‐Herzegovina, Inbreeding Slovakia) and western (Iberian Peninsula, United Kingdom) parts of −0.086 −0.033 −0.091 −0.101 −0.107 −0.178 the distribution (Figure 2), with a contact admixture zone in cen‐ tral Europe (Switzerland). All the additional partitions (K = 3–12) F produced the same results, that is, the separation of Swiss popu‐ lations as a third independent pool and, in some instances, Bosnia‐ Nucleotide diversity per site × 10 −3; All: all sites, including coding regions (CDS), introns, and intergenic sequences; Syn: synonymous sites; Non‐syn: nonsynonymous sites. FuF* 0.57 1.41 1.10 1.11 1.19 0.74 Herzegovina and Slovakia were split into two different genetic pools (data not shown). Neutrality test stats Demographic inference based on stairway plot revealed a de‐ FuD* 1.00 0.54 0.68 1.69 0.90 0.99 crease in effective population size for all populations starting around 400,000–600,000 years ago (Figure 4). Since then, the de‐ mographic decline has been continuous for populations from the 1.09 0.53 1.27 TajD 0.99 0.76 1.19 Iberian Peninsula and United Kingdom (Western Europe yew range), as well as for two of the Swiss populations, reaching current effec‐ tive population sizes (100–400 individuals) around 500–2,000 years 0.308 ago (Figure 4). The remaining Swiss population (WALL) and those 0.312 0.318 0.316 0.313 0.310 πa/πs from the Eastern Europe yew range (Bosnia‐Herzegovina, Slovakia), however, seem to have maintained high (7,000–20,000 individu‐ als) and stable effective population sizes during the last 10,000– Non‐syn 20,000 years (Figure 4). A higher demographic decline for the former 2.54 2.53 2.63 2.63 2.71 2.76 group of populations was also supported by higher positive values of Nucleotide diversitya (π) Tajima's D, and Fu & Li's D* and F*, indicating stronger demographic bottlenecks (Table 1). 8.00 8.40 8.68 8.49 7.99 8.97 Syn 3.3 | Adaptive variation at SNP, gene, and pathway level 4.30 4.37 3.97 3.95 4.16 4.14 All We detected a similar number of significant associations between SNPs and climate with bayescenv and baypass (676 vs. 777), but most 80.38 82.80 78.65 75.27 65.61 of them involved single SNPs per gene. The number of significant 100 Polymorphic sites % associations with at least two different SNPs on the same gene was lower, involving 212/293 SNPs located on 72/111 genes for bayes- cenv and baypass , respectively. However, only 11 SNPs located on 7 Number 20,680 21,302 16,879 20,233 25,726 19,363 genes had significant associations with the same climatic variables using both methods. At gene level, we found 60 genes with a sig‐ nificantly higher proportion of SNPs associated with climate than expected by chance. Combining these results with those obtained 1,472,324 1,472,125 1,470,875 1,472,055 1,473,550 1,471,260 at the SNP level (see Materials and Methods), a total of 11 top can‐ L (bp) didate genes for climate adaptation were retained (Table 2). These potentially adaptive genes were involved in plant growth, develop‐ L, length of sequences; N, sample size. ment, endogenous biological rhythms, and several stress responses (Table 2). 120 40 30 30 10 10 indicate unrelated individuals N The distribution of allelic frequencies in relation to climate of significant SNPs from the 11 top candidate genes confirmed that associations were not due to population structure or demographic Bosnia‐Herzegovina Geographic region Iberian Peninsula processes affecting particular geographic regions (Figure 5 and United Kingdom Figure S6). For instance, for the putative pseudo‐response regulator Switzerland 7 gene (PRR7), a significant higher frequency (F1,10 = 19.13; p < 0.01) Slovakia of the minor allele of SNP5227 was found in those populations All with annual mean precipitation below 1,000 mm, irrespectively of a MAYOL et al. | 151 F I G U R E 4 Demographic history of Taxus baccata populations inferred with stairway plot (Liu & Fu, 2015). Colored lines show the changes in median effective population size for the last one million years, and dotted lines represent the 95% confidence intervals (2.5% and 97.5% percentiles) based on 200 bootstraps. Populations are colored according to their geographic origin: orange (Iberian Peninsula), green (Switzerland), dark pink (United Kingdom), and lilac (eastern populations: Slovakia and Bosnia‐Herzegovina). Colors are the same as in Figure S6. See Table S2 for details on populations their geographic location (Figure 5). The minor allele frequency of Only 390 out of 1,128 genes could be mapped on existing KEGG SNP4561 from early‐responsive to dehydration stress gene (ERD4) in‐ pathway maps for plants. Mapped genes were distributed on 60 dis‐ creased with precipitation of the driest month, being significantly tinct pathways, eight of which were plant‐specific (Figure 6). Only 2 higher in those populations with minimum precipitation above of the 390 genes were included among the 11 top candidate genes 60 mm (F1,10 = 8.50; p < 0.05). Minor allele frequencies of SNP19480 reported above: The putative pseudo‐response regulator 7 gene was (acetyl‐CoA carboxylase 1) and SNP346 (senescence‐associated E3 part of the circadian rhythm‐plant pathway, and the putative acetyl‐ ubiquitin ligase 1) increased significantly in those populations with CoA carboxylase 1 gene was shared between four different path‐ annual mean temperature above 8ºC (F1,10 = 3.81; p = 0.08) and min‐ ways, that is, fatty acid biosynthesis, pyruvate metabolism, propanoate imum temperature of the coldest month above −3°C (F1,10 = 62.24; metabolism, and AMPK signaling pathway. Interestingly, although the p < 0.001), respectively (Figure 5). remaining 388 genes that could be mapped on the pathways showed Similar numbers of environmental correlations were found for no signs of selection when analyzed individually at the SNP or the temperature and precipitation variables (six vs. five, respectively). gene level, we were able to detect a total of seven pathways that Regarding the climatic period, most of the significant correlations de‐ were globally enriched for signals of selection after removing gene tected involved current climate, although a few were also detected redundancy (Q‐value < 0.15; Figure 6). Five pathways (circadian for the last glacial maximum and the last interglacial ones (Table 2). rhythm‐plant, mRNA surveillance, phagosome, spliceosome, and Wnt 152 | TA B L E 2 Candidate genes putatively under selection detected for Taxus baccata. “Gene level” refers to the number of SNPs detected with Yeaman's et al. (2016) approach (see Materials and Methods for details) SNP level Gene Code ConGenIE ID (A. thali‐ (T. baccata) Length (bp) Nº SNPs Description Abbr. ana best BLAST) E‐value Climatic drivers (Period) Gene level BayPass BayeScEnv Tbac_18055 2,773 30 Pseudo‐response regula‐ PRR7 MA_124244g0020 0.0 Annual mean precipitation 15 2 – tor 7 (AT5G02810.1) (PRE) Tbac_19823 1,650 16 P‐loop containing nu‐ DNK MA_892200g0010 2.0E‐170 Annual mean precipitation 7 2 – cleoside triphosphate (AT1G72040.1) (PRE) hydrolases superfamily protein Tbac_129399 843 6 Early‐responsive to de‐ ERD4 MA_99251g0010 4.3E‐59 Precipitation of the driest 5 (PRE) 4 (PRE) 4 (PRE)a hydration stress protein (AT1G30360.1) month (PRE, LIG) 5 (LIG) 4 (LIG) 4 (LIG)a (ERD4) Tbac_71106 951 12 TATA‐binding protein‐as‐ BTAF1 MA_79681g0010 1.8E‐151 Precipitation of the driest 7 (PRE) 2 (PRE) ‐ sociated factor BTAF1 (AT3G54280.1) month (PRE, LIG) 7 (LIG) 3 (LIG) Tbac_27725 856 19 Plant viral‐response family DUF716 MA_172105g0010 9.2E‐111 Annual mean precipitation 9 3 – protein (DUF716) (AT5G13890.3) (PRE) Precipitation of the wettest 10 2 – month (PRE) Tbac_70534 4,869 33 Acetyl‐CoA carboxylase 1 ACC1 MA_38431g0020 1.5E‐83 Annual mean temperature 11 4 – (AT1G36160.2) (PRE) Tbac_19996 1,978 42 Putrescine‐binding peri‐ ENF2 MA_10434660g0010 0.0 Maximum temperature of the 18 4 8 plasmic protein‐related (AT1G31410.1) warmest month (PRE) Tbac_72122 947 16 AAA‐type ATPase family 3A‐ATP MA_8790100g0010 5.1E‐59 Maximum temperature of the 7 2 – protein (AT2G45500.2) warmest month (LGM) Tbac_101255 985 13 Alternative oxidase 1A AOX1A MA_10430050g0010 5.6E‐141 Annual mean temperature 7 3 1 (AT3G22370.1) (LGM) 7 4 1 Minimum temperature of the coldest month (LGM) Tbac_101139 1,306 22 Senescence‐associated E3 SAUL1 MA_113300g0020 7.6E‐129 Annual mean temperature 8 6 – ubiquitin ligase 1 (AT1G20780.1) (LIG) Minimum temperature of the 11 8 – coldest month (PRE) Tbac_70556 6,025 87 Modifier of SNC1 MOS1 MA_958g0010 0.0 Minimum temperature of the 38 2 – (AT4G24680.1) coldest month (LGM) Note: LGM, Last Glacial Maximum; LIG, Last Interglacial, PRE, present. a Q‐values < 0.07 for precipitation of the driest month (PRE, LIG); bayescenv also supported a role of annual mean temperature during PRE and LIG (Q‐values < 0.05). MAYOL et al. MAYOL et al. | 153 F I G U R E 5 Geographic distribution of minor allele frequency for single SNPs of four top candidate genes and their relation to climate. ERD4: early‐responsive to dehydration stress protein (SNP4561, G/A SNP, A allele plotted); PRR7: pseudo‐response regulator 7 (SNP5227, T/C SNP, C allele plotted); ACC1: acetyl‐CoA carboxylase 1 (SNP19480, G/C SNP, G allele plotted); SAUL1: senescence‐associated E3 ubiquitin ligase 1 (SNP346, T/C SNP, C allele plotted). The size of the circle is proportional to allele frequency. Precipitation and temperature variables are shown in blue and red, respectively. Darker colors indicate higher values of climatic variables signaling) showed significantly lower values of πa/πs relative to the at the single‐locus (SNP) and multilocus (gene, pathway) levels. Our total number of pathways analyzed, indicating constrained adaptive results also show that polygenic approaches can provide useful and evolution (Table 3). This was because many genes (>25%) in these simple tools for detecting signatures of “soft” selection (i.e., acting pathways had no variation for nonsynonymous mutations, that is, on multiple genes at the same time). Using this polygenic approach, were highly conserved (Table S5). The Oxidative phosphorylation we were able to identify one biological pathway as potentially under pathway had significantly lower global values of π, pointing also to the effect of diversifying selection, when individual approaches such conserved genes, and flavonoid biosynthesis exhibited significantly as those focused on single loci failed to detect the signature of selec‐ higher values of π and Fu & Li's D*, a signature of widespread diver‐ tion for any of the genes in this pathway. Finally, our study provides sifying selection in this pathway (Table 3; Table S5). None of these some evidence supporting the contribution of past climate to the seven gene pathways contained genes with extreme values causing potentially adaptive genetic variation observed. Below we discuss highly deviant SUMSTAT scores, but they had significantly lower or these findings within the context of the demography of the species higher values because all genes in the pathway globally shifted to while considering their potential limitations and caveats. small or large values (Table S5). 4.1 | Evidence of selection in response to past and recent climate despite demographic decline 4 | D I S CU S S I O N Using extensive genetic resources (1,128 genes and 25,726 SNPs), Our study demonstrates the potential advantage of using multiple we confirmed previous demographic studies on English yew that complementary approaches for detecting selection to overcome supported the existence of two main Eastern and Western genetic limitations inherent in testing them individually. For instance, the pools across the whole species range (Mayol et al., 2015). In addi‐ clinal pattern found for phenological traits (i.e., male strobili mat‐ tion, we found a strong signal of demographic decline that was highly uration, proportion of late growth) suggests the existence of local consistent across populations, suggesting that common factors adaptations at the Iberian scale, and this assumption is reinforced may underlie demographic processes in this species and confirm‐ by the detection of similar signatures at the European scale, both ing on a wider scale previous reports from the Iberian populations 154 | MAYOL et al. F I G U R E 6 Diagram of the 60 pathways tested for signals of polygenic selection (* indicates pathways that were plant‐specific). The size of the nodes (pathways) is proportional to the number of genes (minimum size = 5, maximum size = 19). Edges represent shared genes between pathways, and edge width is proportional to the number of shared genes (minimum width = 1, maximum width = 7). Candidate pathways for polygenic adaptation are indicated in red. Colored circles indicate significant lower or higher global statistics (Q‐values < 0.15; see text and Table 3 for further details). Different colors refer to different statistics: orange to πa/πs (ratio of nonsynonymous to synonymous nucleotide diversity); purple to nucleotide diversity (π); green to both π and Fu & Li's D* TA B L E 3 Candidate pathways for polygenic adaptation in Taxus baccata after removing overlapping genes (“pruning,” see Daub et al., 2013) Pathway size before/after Q‐value (only Gene pathway pruning Statistics tested p‐value values < 0.15) circadian rhythm‐plant 6/5 πa/πs 0.005619 0.083 mRNA surveillance 7/7 πa/πs 0.003249 0.082 phagosome 8/6 πa/πs 0.004259 0.082 spliceosome 10/10 πa/πs 0.003789 0.082 Wnt signaling 10/10 πa/πs 0.000009 0.003 Oxidative phosphorylation 9/9 π 0.003099 0.073 flavonoid biosynthesis 10/10 π 0.003259 0.073 flavonoid biosynthesis 10/5 FuD* 0.000239 0.011 FuD*: Fu & Li's D* (Fu & Li, 1993). (Burgarella et al., 2012). Although the inferred start of population Beaulieu et al., 2001; Koutsodendris et al., 2010). We found, how‐ decline (~400,000–600,000 years ago) was somewhat earlier than ever, that one population from Switzerland and the two from the that reported by Burgarella et al. (2012), who placed it at ~100,000– eastern range (Bosnia‐Herzegovina, Slovakia) were characterized by 300,000 years ago, both studies show the same trend for the spe‐ stable and larger population sizes, which is in agreement with the cies as a whole, which suggests that habitat suitability for English pattern of decreasing genetic diversity from east to west reported yew has been worsening for the last several millennia, a pattern for English yew, and could be explained by the milder climatic condi‐ that is also evident in most European pollen records since the last tions experienced by easternmost populations since the long past interglacial period (c. 115,000–130,000 year BP; Turner, 2000; de (Mayol et al., 2015). MAYOL et al. | 155 Interestingly, despite the small current effective population important differences exist between angiosperms and gymno‐ size inferred for most populations, we found evidence for poten‐ sperms clock systems (Gyllenstrand et al., 2014). tially adaptive variation at both the phenotypic and molecular level. Interestingly, a significant association with latitude has been re‐ Following theoretical considerations, Kawecki and Ebert (2004) ported for candidate SNPs at the PRR7 gene in European popula‐ suggested that convincing evidence for local adaptation requires tions of Norway spruce (Chen, Källman, et al., 2012), in accordance testing for population × habitat interaction in reciprocal transplant with the pattern detected for English yew and annual mean precip‐ experiments. However, even a single common garden experiment itation at the European scale. Although most studies investigating can help to reveal genetic differentiation among widely distributed, the adaptive response of circadian genes have focused on the role of divergent populations, providing indirect evidence of past adapta‐ temperature and photoperiod, recent studies in Arabidopsis indicate tion (Cheplick, 2015). In our study, the differential growth (a proxy that daily humidity oscillations are also involved in the regulation of for fitness in forest trees) of the Iberian yew populations was as‐ specific circadian outputs besides those co‐regulated with the light– sociated with the climatic distance between the site of origin and dark cycle (Mwimba et al., 2018). Moreover, the PRR7 gene seems to the conditions experienced in the common environment, suggest‐ be directly involved in the regulation of both cold and drought stress ing local adaptation. Similarly, the among‐population variability in responses of plants in coordination with an ABA‐dependent mech‐ shoot elongation and male strobili maturation, both linked to tem‐ anism (Liu, Carlsson, Takeuchi, Newton, & Farré, 2013). This could perature clines, also pointed to the selective effects of climate on explain the apparent contradiction between the results obtained in phenological patterns. Compared to more temperate coastal pop‐ the common garden, that probably reflect the overall response of ulations, those from continental environments (colder and subject the circadian clock genes to temperature, and the variation found to more pronounced temperature variations) showed both faster for the PRR7 gene at the European scale in response to precipitation. growth and earlier maturation of male strobili at the trial site, which Apart from PRR7, several other candidate genes relevant for for‐ again suggests that natural selection has favored genotypes opti‐ est tree adaptation were detected in our study (Table 2). Some of mized to the short time frame available for successful growth and these genes are rapidly activated in response to drought stress (early‐ reproduction. Of course, part of the observed variability among responsive to dehydration stress (ERD4): Alves & Fietto, 2013), cold genets and populations could also be due to strong maternal clone acclimation (acetyl‐CoA carboxylase 1 (ACC1): Amid, Lytovchenko, effects (C‐effects, see Burdon & Shelbourne, 1974). These transgen‐ Fernie, Warren, & Thorlby, 2012; Alternative oxidase 1A (AOX1A): erational clonal effects have been widely reported (see review in Fiorani, Umbach, & Siedow, 2005; Putrescine‐binding periplasmic Schwaegerle, Mcintyre, & Swingley, 2000). However, experimental protein‐related (ENF2): Cuevas et al., 2008), and pathogen effectors treatments in plants suggest that the magnitude of C‐effects usually (Modifier of SNC1 (MOS1): Li, Tessaro, Li, & Zhang, 2010), while oth‐ vanishes with time, particularly for late developing traits (e.g., mass ers are central components of the regulatory network controlling or growth; Libby & Jund, 1962; Schwaegerle et al., 2000). Therefore, leaf senescence in plants (senescence‐associated E3 ubiquitin ligase while C‐effects cannot be fully ruled out, their importance is ex‐ 1 (SAUL1): Lee et al., 2016). However, despite the candidate genes pected to be relatively low, since measures of plant performance reported in this study provide a first set of genes to further explore were obtained from plants growing on a common environment for English yew's adaptive responses at the molecular level, there is a 10–15 years. Many studies on tree species have also reported clinal need to confirm our results using a higher number of samples, as well patterns of phenotypic differentiation, in particular for phenolog‐ as to provide unequivocal links between allelic variation and fitness ical and fitness‐related traits, suggesting that tree populations are before concluding that the patterns found are truly adaptive. Both locally adapted over large parts of their native ranges (see review theoretical and empirical studies suggest that accurate estimates in Savolainen, Pyhäjärvi, & Knürr, 2007; Alberto et al., 2013). Our of genetic differentiation can be obtained with a large number of findings for English yew are similar to these studies but, unlike many SNPs (≥1,500) even at small sample sizes (Nazareno, Bemmels, Dick, widespread “model” tree species, adaptation occurs in the face & Lohmann, 2017; Willing, Dreyer, & Van Oosterhout, 2012) but, in of limited effective population sizes and restricted gene flow (see some cases, detection of regions under selection may be affected by Dubreuil et al., 2010). nonrandom processes such as population structuring. Detecting the Analyses at the molecular level also suggested local adaptation same patterns in a larger number of samples could help to circum‐ for phenological traits, since significant signatures of selection were vent these problems and to increase the confidence in our results. found at the European scale for both the circadian rhythm pathway Although most of the variation found for these candidate genes (discussed further in the next section) and the putative pseudo‐re‐ was correlated to present climate, we also were able to detect sig‐ sponse regulator 7 gene (PRR7). In A. thaliana, PRR7 is an essential nificant associations with past climate for a few genes (Table 2). It component of the circadian clock that is directly associated with de‐ must be noted, however, that our approximation to detect the po‐ velopmental transitions such as flowering (Nakamichi et al., 2007), tential contribution of past climate is not free of uncertainty, since and putative orthologous genes have been identified in Picea abies it is based on the assumption that current distribution of English (L.) Karst. (Gyllenstrand et al., 2014). In this conifer, PRR7 shows yew has been similar across the different Quaternary periods. On partly conserved protein domains and displays diurnal rhythms the other hand, the severe bottleneck that the populations have of expression similar to those observed in Arabidopsis, although undergone since the long past probably has contributed to the loss 156 | MAYOL et al. of variants, making more difficult the detection of signatures of se‐ the whole pathway, the PRR7 gene seems to be a good candidate lection in response to past climate. Despite these potential caveats, for adaptive variation in response to climate in English yew, as dis‐ there is some evidence suggesting that some of the associations cussed in the previous section. These results are highly consistent found for allelic variation and past climates could be adaptive. For with those found in Picea abies, where both purifying and diversify‐ instance, none of the associations found for the current climate were ing selection have been invoked to explain the diversity found for also found for the LGM. The few genes (three) whose variation was genes in the photoperiodic (circadian rhythm) pathway (Chen et al., found to be linked to more than one climatic period involved both 2016). In this species, upstream light receptor genes, such as cryp‐ interglacials (PRE, LIG), with two of them also associated with the tochrome‐1, phytocromes, and Zeitlupe, have about half the nucleo‐ same environmental variable (i.e., precipitation of the driest month; tide diversity than the downstream genes of the circadian clock, like Table 2). This could just reflect a higher climatic similarity between Gigantea and the pseudo‐response regulators 1, 3, and 7 (π = 0.0016 the present day climate and that of the last interglacial, but also vs. π = 0.0031, respectively; Källman et al., 2014). Accordingly, the points to a differential response to glacial and interglacial climates. signals of diversifying selection are mostly identified in downstream Additional evidence is provided by the lack of correlation between genes (Chen, Källman, et al., 2012; Chen et al., 2016; Källman et allelic variation and precipitation variables during LGM (see Table 2), al., 2014). These results are highly consistent with our findings for suggesting that water availability could have acted as a more im‐ T. baccata (Table S5), where the nucleotide diversity for the puta‐ portant selective pressure during interglacials than during glacial tive PRR7 gene (π = 0.0036) was 3‐fold that of the remaining genes periods. In contrast, similar temperature‐mediated responses were (π = 0.0011). In addition, the ratio of nonsynonymous to synony‐ found for cold and warm glacial episodes, since the same number of mous variation for the PRR7 gene was between 5–30 times greater candidate genes involved in thermal adaptation was detected for the than that of the other genes in this pathway, suggesting that func‐ glacial (LGM) and interglacial (PRE, LIG) periods. However, further tional constraints for this particular gene might be much more exploration of the relationship between adaptive variation and past relaxed (Table S5). In fact, results for both P. abies and T. baccata and current climates in other tree species will be needed to fully test are in agreement with the prediction that selective constraint is this hypothesis. progressively relaxed along metabolic pathways because upstream genes might have greater pleiotropic effects than downstream ones, since upstream genes are required for a wider range of end 4.2 | Signature of polygenic adaptation in products (Cork & Purugganan, 2004), as has been proven in other biological pathways metabolic pathways in plants (e.g., carotenoid biosynthesis in carrot, In addition to single‐locus methods, we also applied a gene set en‐ Clotaut, Peltier, Soufflet‐Freslon, Briard, & Geoffriau, 2012; isofla‐ richment test that has been used successfully to detect polygenic vonoid pathway in soybean, Chu, Wang, Cheng, Yang, & Yu, 2014; adaptation in the human genome (Daub et al., 2013, 2017; Foll et but see Olson‐Manning, Lee, Rausher, & Mitchell‐Olds, 2013 for al., 2014). With this approach, we found seven biological pathways glucosinolate pathway in Arabidopsis thaliana). that were globally enriched for signals of selection, being therefore Finally, the flavonoid biosynthesis pathway had significantly good candidates for polygenic adaptation. Five pathways involved in higher values for overall nucleotide diversity (π) and Fu & Li's D*, essential primary functions that are critical for plant survival, namely pointing to divergent selection acting across English yew's geo‐ catabolism (phagosome), genetic information processing (mRNA graphic range. Flavonoids are plant secondary metabolites that are surveillance, spliceosome) or environmental information processing involved in a wide array of biological functions, including protection (circadian rhythm, Wnt signaling), were characterized by highly con‐ against ultraviolet (UV) radiation and phytopathogens, signaling served genes (that is, significantly lower overall values of the ratio of during nodulation, pollen fertility, auxin transport and coloration of nonsynonymous to synonymous nucleotide diversity). In fact, many flowers, fruits and seeds (Falcone Ferreyra, Rius, & Casati, 2012). of the constituent genes of these pathways showed no variation for Interestingly, neither analysis at the SNP nor at the single‐gene level nonsynonymous mutations (Table S5), suggesting that they evolved showed significant signatures of selection for any of the genes in under rigorous functional or structural constraints. An additional this pathway. Thus, our approach of testing for selection at the path‐ pathway, which is involved in crucial energy metabolism processes way level may have allowed us to identify genetic responses to the (Oxidative phosphorylation), showed significant lower overall values local environment as a whole, not only to climate (as in the SNP/ of nucleotide diversity. These results are consistent with the ex‐ gene‐level approaches) and/or better assess signatures of selection pected prevalence of purifying selection for genes and/or genomic for highly polygenic adaptive traits. Several genes with small effect regions of functional relevance, especially for nonsynonymous sites mutations can together have a large impact on a biological pathway (Wright & Andolfatto, 2008), driving a global shift toward low values and, as stress‐responsive genes, those coding for flavonoids are ex‐ of nucleotide diversity and/or high efficacy of selection (that is, low pected to evolve under the distinct selection pressures of each local πa/πs). environment, including day length, light quality, water availability Within this group, however, the circadian rhythm pathway de‐ and temperature, as well as local pathogens. It is thus possible that serves special attention (Figure S7). Despite the prevalent role of the higher overall nucleotide diversity and Fu & Li's D* values reflect purifying selection suggested by the gene enrichment approach for slight changes of allele frequency but at many genes at the same MAYOL et al. | 157 time, which is difficult to identify with conventional single‐locus C O N FL I C T O F I N T E R E S T approaches (Pritchard et al., 2010). This opens a promising line of None declared. research to be applied to other nonmodel forest trees. DATA AC C E S S I B I L I T Y 4.3 | Concluding remarks Data for this study are available at Dryad Digital Repository: https​:// Our rangewide study in English yew using a large number of SNP doi.org/10.5061/dryad.78j501j (Mayol et al., 2019). markers allowed us to confirm an overall demographic history characterized by two gene pools in Western and Eastern Europe, with an admixture zone in Central Europe, and population decline ORCID in most populations, extending initial evidence based on nuSSRs in Maria Mayol https://orcid.org/0000-0001-8407-9083 Burgarella et al. (2012). Interestingly, despite demographic decline, we found evidence of local adaptation to climate both for phenotypic Miquel Riba https://orcid.org/0000-0003-2164-9880 traits in a common garden and SNPs from several candidate genes, Stephen Cavers https://orcid.org/0000-0003-2139-9236 probably because levels of nucleotide diversity (π = 0.00437) are still Delphine Grivet https://orcid.org/0000-0001-8168-4456 large in this long‐lived tree. Several new candidate genes for positive Lucie Vincenot https://orcid.org/0000-0002-6481-0139 and negative selection in forest trees were identified or confirmed in this study, and an approach to detect potential polygenic selection in Federica Cattonaro https://orcid.org/0000-0001-8819-7458 biological pathways was successfully tested. In particular, we identi‐ Giovanni G. Vendramin https://orcid.org/0000-0001-9921-7872 fied the flavonoid biosynthesis pathway as a general stress‐response Santiago C. González‐Martínez https://orcid. pathway that deserves further attention, including genomic studies org/0000-0002-4534-3766 in other plant systems. Our study also suggests that the integration of complementary information coming from different sources could be useful to detect selection in nonmodel species for which com‐ REFERENCES mon gardens and genomic resources are limited. Such an approach Ahrens, C. W., Rymer, P. D., Stow, A., Bragg, J., Dillon, S., Umbers, K. D. L., can also help the field to advance beyond correlational “black‐box” & Dudaniec, R. Y. (2018). The search for loci under selection: Trends, methods (i.e., “reverse ecology”) to systems biology‐driven ones, biases and progress. Molecular Ecology, 27, 1342–1356. https​://doi. which are more informative in terms of drivers and mechanisms of org/10.1111/mec.14549​ adaptation. Finally, our study also contributed to the emerging view Ahrens, J. H., & Dieter, U. (1985). Sequential random sampling. ACM Transactions on Mathematical Software, 11, 157–169. https​://doi. that while adaptation to current climates explains part of the con‐ org/10.1145/214392.214402 temporary distribution of standing genetic variation, it is also essen‐ Alberto, F. J., Aitken, S. N., Alía, R., González‐Martínez, S. C., Hänninen, tial to consider adaptation to past climates, especially for long‐lived H., Kremer, A., … Savolainen, O. (2013). Potential for evolutionary re‐ forest trees. Understanding the drivers and mechanisms of adapta‐ sponses to climate change – evidence from tree populations. Global Change Biology, 19, 1645–1661. https​://doi.org/10.1111/gcb.12181​ tion is particularly relevant for highly threatened species such as the Alía, R., Chambel, R., Notivol, E., Climent, J., & González‐Martínez, S. C. English yew, which is protected by environmental laws in several (2014). Environment‐dependent microevolution in a Mediterranean European countries, and may help to develop sound conservation pine (Pinus pinaster Aiton). BMC Evolutionary Biology, 14, 200. https​:// programs and reforestation schemes. doi.org/10.1186/s12862-014-0200-5 Altschul, S. F., Gish, W., Miller, W., Myers, E. W., & Lipman, D. J. (1990). Basic local alignment search tool. Journal of Molecular Biology, 215, AC K N OW L E D G E M E N T S 403–410. https​://doi.org/10.1016/S0022-2836(05)80360-2 Alves, M. S., & Fietto, L. G. (2013). Functional diversity of early respon‐ We kindly thank D. Ballian and L. Paule for providing yew samples. We sive to dehydration (ERD) genes in soybean. In J. E. Board (Ed.), A comprehensive survey of international soybean research – Genetics, are also indebted to E. Berganzo for her support with clonal bank meas‐ Physiology, Agronomy and Nitrogen Relationships (pp. 475–487). urements, C. García‐Barriga for laboratory assistance, and I. Regalado Rijeka, Croatia: IntechOpen. ISBN 978-953-51-0876-4. http://dx.doi. and C. Batlles for help in sample collection. We thank the reviewers org/10.5772/45867​ for their helpful comments and for the suggestion of a new title. This Amid, A., Lytovchenko, A., Fernie, A. R., Warren, G., & Thorlby, G. J. work was supported by projects TAXUS (CGL2007‐63107⁄BOS) and (2012). The sensitive to freezing3 mutation of Arabidopsis thaliana is a cold‐sensitive allele of homomeric acetyl‐CoA carboxylase that ADAPCON (CGL2011‐30182‐C02‐01/02) from the Spanish Ministry results in cold‐induced cuticle deficiencies. Journal of Experimental of Economy and Competitiveness. This publication is also part of a Botany, 63, 5289–5299. https​://doi.org/10.1093/jxb/ers191 project funded by the European Union's Horizon 2020 research and Barrett, R. D. H., & Hoekstra, H. E. (2011). Molecular spandrels: Tests of innovation program under grant agreement No. 676876 (GenTree). adaptation at the genetic level. Nature Reviews Genetics, 12, 767–780. https​://doi.org/10.1038/nrg3015 Special thanks to S. Iglesias and the Spanish Ministry of Agriculture, Bates, D., Maechler, M., Bolker, B., & Walker, S. (2015). Fitting linear Fisheries, and Food for granting permit to work in the Valsaín Clonal mixed‐effects models using lme4. Journal of Statistical Software, 67, Bank and all related logistic support. 1–48. 158 | MAYOL et al. Berg, J. J., & Coop, G. (2014). A population genetic signal of polygenic primate lineages ancestral to humans. Molecular Biology & Evolution, adaptation. PLoS Genetics, 10, e1004412. https​://doi.org/10.1371/ 34, 1391–1402. https​://doi.org/10.1093/molbe​v/msx083 journ​al.pgen.1004412 de Beaulieu, J. L., Andrieu‐Ponel, V., Reille, M., Grüger, E., Tzedakis, C., Burdon, R. D., & Shelbourne, C. J. A. (1974). The use of vegetative propa‐ & Svobodova, H. (2001). An attempt at correlation between the gules for obtaining genetic information. New Zealand Journal of Forest Velay pollen sequence and the Middle Pleistocene stratigraphy from Science, 4, 418–425. Central Europe. Quaternary Science Reviews, 20, 1593–1602. https​:// Burgarella, C., Navascués, M., Zabal‐Aguirre, M., Berganzo, E., doi.org/10.1016/S0277-3791(01)00027-0 Riba, M., Mayol, M., … González‐Martínez, S. C. (2012). Recent Del Fabbro, C., Scalabrin, S., Morgante, M., & Giorgi, F. M. (2013). An population decline and selection shape diversity of taxol‐re‐ extensive evaluation of read trimming effects on Illumina NGS lated genes. Molecular Ecology, 21, 3006–3021. https​://doi. data analysis. PLoS ONE, 8, e85024. https​://doi.org/10.1371/journ​ org/10.1111/j.1365-294X.2012.05532.x al.pone.0085024 Chen, J., Glemin, S., & Lascoux, M. (2017). Genetic diversity and the effi‐ Dubreuil, M., Riba, M., González‐Martíınez, S. C., Vendramin, G. G., cacy of purifying selection across plant and animal species. Molecular Sebastiani, F., & Mayol, M. (2010). Genetic effects of chronic hab‐ Biology & Evolution, 34, 1417–1428. https​://doi.org/10.1093/molbe​v/ itat fragmentation revisited: Strong genetic structure in a temper‐ msx088 ate tree, Taxus baccata (Taxaceae), with great dispersal capability. Chen, J., Källman, T., Ma, X.‐F., Gyllenstrand, N., Zaina, G., Morgante, American Journal of Botany, 97, 303–310. M., … Lascoux, M. (2012). Disentangling the roles of history and Epps, C. W., & Keyghobadi, N. (2015). Landscape genetics in a chang‐ local selection in shaping clinal variation of allele frequencies and ing world: Disentangling historical and contemporary influences gene expression in Norway Spruce (Picea abies). Genetics, 191, and inferring change. Molecular Ecology, 24, 6021–6040. https​://doi. 865–881. org/10.1111/mec.13454​ Chen, J., Källman, T., Ma, X.‐F., Zaina, G., Morgante, M., & Lascoux, M. Falcone Ferreyra, M. L., Rius, S. P., & Casati, P. (2012). Flavonoids: (2016). Identifying genetic signatures of natural selection using pooled Biosynthesis, biological functions, and biotechnological applications. population sequencing in Picea abies. G3: Genes|genomes|genetics, 6, Frontiers in Plant Science, 3, 222. 1979–1989. Fiorani, F., Umbach, A. L., & Siedow, J. N. (2005). The alternative oxidase Chen, J., Uebbing, S., Gyllenstrand, N., Lagercrantz, U., Lascoux, M., & of plant mitochondria is involved in the acclimation of shoot growth Källman, T. (2012). Sequencing of the needle transcriptome from at low temperature. A study of Arabidopsis AOX1a transgenic plants. Norway spruce (Picea abies Karst L.) reveals lower substitution rates, Plant Physiology, 139, 1795–1805. but similar selective constraints in gymnosperms and angiosperms. Foll, M., Gaggiotti, O. E., Daub, J. T., Vatsiou, A., & Excoffier, L. (2014). BMC Genomics, 13, 589. Widespread signals of convergent adaptation to high altitude in Asia Cheplick, G. P. (2015). Approaches to plant evolutionary ecology. (pp. 46). and America. The American Journal of Human Genetics, 95, 394–407. New York, NY: Oxford University Press. https​://doi.org/10.1016/j.ajhg.2014.09.002 Chu, S., Wang, J., Cheng, H., Yang, Q., & Yu, D. (2014). Evolutionary study Fu, Y. X., & Li, W. H. (1993). Statistical tests of neutrality of mutations. of the isoflavonoid pathway based on multiple copies analysis in soy‐ Genetics, 133, 693–709. bean. BMC Genetics, 15, 76. https​://doi.org/10.1186/1471-2156-15-76 Gautier, M. (2015). Genome‐wide scan for adaptive divergence and as‐ Chybicki, I. J., Oleksa, A., & Kowalkowska, K. (2012). Variable rates sociation with population‐specific covariates. Genetics, 201, 1555– of random genetic drift in protected populations of English yew: 1579. https​://doi.org/10.1534/genet​ics.115.181453 Implications for gene pool conservation. Conservation Genetics, 13, González‐Martínez, S. C., Dubreuil, M., Riba, M., Vendramin, G. G., 899–911. https​://doi.org/10.1007/s10592-012-0339-9 Sebastiani, F., & Mayol, M. (2010). Spatial genetic structure of Taxus Clotault, J., Peltier, D., Soufflet‐Freslon, V., Briard, M., & Geoffriau, E. baccata L. in the western Mediterranean Basin: Past and present (2012). Differential selection on carotenoid biosynthesis genes as a limits to gene movement over a broad geographic scale. Molecular function of gene position in the metabolic pathway: A study on the Phylogenetics and Evolution, 55, 805–815. carrot and dicots. PLoS ONE, 7, e38724. Gugger, P. F., Ikegami, M., & Sork, V. L. (2013). Influence of late Quaternary Cork, J. M., & Purugganan, M. D. (2004). The evolution of molecular ge‐ climate change on present patterns of genetic variation in valley oak, netic pathways and networks. BioEssays, 26, 479–484. https​://doi. Quercus lobata Née. Molecular Ecology, 22, 3598–3612. org/10.1002/bies.20026​ Gyllenstrand, N., Karlgren, A., Clapham, D., Holm, K., Hall, A., Gould, P. Csillérry, K., Lalagüe, H., Vendramin, G. G., González‐Martínez, S. C., D., … Lagercrantz, U. (2014). No time for spruce: Rapid dampening of Fady, B., & Oddou‐Muratorio, S. (2014). Detecting short spatial scale circadian rhythms in Picea abies (L. Karst). Plant & Cell Physiology, 55, local adaptation and epistatic selection in climate‐related candidate 535–550. https​://doi.org/10.1093/pcp/pct199 genes in European beech (Fagus sylvatica) populations. Molecular Hao, D. C., Ge, G., Xiao, P., Zhang, Y. Y., & Yang, L. (2011). The first insight Ecology, 23, 4696–4708. into the tissue specific Taxus transcriptome via Illumina second gen‐ Cuevas, J. C., López‐Cobollo, R., Alcázar, R., Zarza, X., Koncz, C., Altabella, eration sequencing. PLoS ONE, 6, e21220. T., … Ferrando, A. (2008). Putrescine is involved in Arabidopsis freez‐ Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., & Jarvis, A. (2005). ing tolerance and cold acclimation by regulating abscisic acid levels Very high resolution interpolated climate surfaces for global land in response to low temperature. Plant Physiology, 148, 1094–1105. areas. International Journal of Climatology, 25, 1965–1978. https​:// https​://doi.org/10.1104/pp.108.122945 doi.org/10.1002/joc.1276 Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, Holliday, J. A., Aitken, S. N., Cooke, J. E. K., Fady, B., González‐Martínez, M. A., … 1000 Genomes Project Analysis Group (2011). The variant S. C., Heuertz, M., … Plomion, C. (2017). Advances in ecological ge‐ call format and VCFtools. Bioinformatics, 27, 2156–2158. https​://doi. nomics in forest trees and applications to genetic resources con‐ org/10.1093/bioin​forma​tics/btr330 servation and breeding. Molecular Ecology, 26, 706–717. https​://doi. Daub, J. T., Hofer, T., Cutivet, E., Dupanloup, I., Quintana‐Murci, L. L., org/10.1111/mec.13963​ Robinson‐Rechavi, M., & Excoffier, L. (2013). Evidence for polygenic Källman, T., De Mita, S., Larsson, H., Gyllenstrand, N., Heuertz, M., adaptation to pathogens in the human genome. Molecular Biology & Parducci, L., … Lascoux, M. (2014). Patterns of nucleotide diver‐ Evolution, 30, 1544–1558. https​://doi.org/10.1093/molbe​v/mst080 sity at photoperiod related genes in Norway spruce [Picea abies Daub, J. T., Moretti, S., Davydov, I. I., Excoffier, L., & Robinson‐Rechavi, (L.) Karst.]. PLoS ONE, 9, e953069. https​://doi.org/10.1371/journ​ M. (2017). Detection of pathways affected by positive selection in al.pone.0095306 MAYOL et al. | 159 Kawecki, T. J., & Ebert, D. (2004). Conceptual issues in local from an Amazonian plant species. Molecular Ecology Resources, 17, adaptation. Ecology Letters, 7, 1225–1241. https​://doi. 1136–1147. https​://doi.org/10.1111/1755-0998.12654​ org/10.1111/j.1461-0248.2004.00684.x Ninyerola, M., Pons, X., & Roure, J. M. (2005). Atlas Climático Digital de Koutsodendris, A., Müller, U. C., Pross, J., Brauer, A., Kotthoff, U., & la Península Ibérica. Metodología y aplicaciones en bioclimatología y Lotter, A. F. (2010). Vegetation dynamics and climate variability geobotánica. Bellaterra: Universidad Autónoma De Barcelona. ISBN during the Holsteinian interglacial based on a pollen record from 932860‐8‐7. Dethlingen (northern Germany). Quaternary Science Reviews, 29, Nystedt, B., Street, N. R., Wetterbom, A., Zuccolo, A., Lin, Y.‐C., Scofield, 3298–3307. https​://doi.org/10.1016/j.quasc​irev.2010.07.024 D. G., … Jansson, S. (2013). The Norway spruce genome sequence Lafontaine, G., Napier, J. D., Petit, R. J., & Hu, F. S. (2018). Invoking adap‐ and conifer genome evolution. Nature, 497, 579–584. https​://doi. tation to decipher the genetic legacy of past climate change. Ecology, org/10.1038/natur​e12211 99, 1530–1546. https​://doi.org/10.1002/ecy.2382 Olson‐Manning, C. F., Lee, C.‐R., Rausher, M. D., & Mitchell‐Olds, T. Lee, E.‐K., Jin, Y.‐W., Park, J. H., Yoo, Y. M., Hong, S. M., Amir, R., … Loake, (2013). Evolution of flux control in the glucosinolate pathway in G. J. (2010). Cultured cambial meristematic cells as a source of plant Arabidopsis thaliana. Molecular Biology and Evolution, 30, 14–23. https​ natural products. Nature Biotechnoly, 28, 1213–1217. https​://doi. ://doi.org/10.1093/molbe​v/mss204 org/10.1038/nbt.1693 Olsson, S., Pinosio, S., González‐Martínez, S. C., Abascal, F., Mayol, M., Lee, I. H., Lee, I. C., Kim, J., Kim, J. H., Chung, E. H., Kim, H. J., … Lim, P. O. Grivet, D., & Vendramin, G. G. (2018). De novo assembly of English (2016). NORE1/SAUL1 integrates temperature‐dependent defense yew (Taxus baccata) transcriptome and its applications for intra‐ and programs involving SGT1b and PAD4 pathways and leaf senescence inter‐specific analyses. Plant Molecular Biology, 97, 337–345. https​:// in Arabidopsis. Physiologia Plantarum, 158, 180–199. doi.org/10.1007/s11103-018-0742-9 Length, R. (2018). emmeans: Estimated marginal means, aka least‐squares Pritchard, J. K., Pickrell, J. K., & Coop, G. (2010). The genetics of human means. R Package Version 1.2. https​://CRAN.R–proje​c t.org/packa​ adaptation: Hard sweeps, soft sweeps, and polygenic adapta‐ ge=emmeans tion. Current Biology, 20, R208–R215. https​://doi.org/10.1016/j. Li, Y. F., Costello, J. C., Holloway, A. K., & Hahn, M. W. (2008). “Reverse cub.2009.11.055 ecology” and the power of population genomics. Evolution, 62, 2984– R Core Team (2018). R: A language and environment for statistical comput‐ 2994. https​://doi.org/10.1111/j.1558-5646.2008.00486.x ing. Vienna, Austria: R Foundation for Statistical Computing. Li, Y., Tessaro, M. J., Li, X., & Zhang, Y. (2010). Regulation of the expres‐ Raj, A., Stephens, M., & Pritchard, J. K. (2014). fastSTRUCTURE: sion of plant resistance gene SNC1 by a protein with a conserved Variational inference of population structure in large SNP data BAT2 domain. Plant Physiology, 153, 1425–1434. sets. Genetics, 197, 573–589. https​://doi.org/10.1534/genet​ Libby, W. J., & Jund, E. (1962). Variance associated with cloning. Heredity, ics.114.164350 17, 533–540. https​://doi.org/10.1038/hdy.1962.55 Rutter, M. A., & Fenster, C. B. (2007). Testing for adaptation to climate in Liu, T., Carlsson, J., Takeuchi, T., Newton, L., & Farré, E. M. (2013). Arabidopsis thaliana: A calibrated common garden approach. Annals Direct regulation of abiotic responses by the Arabidopsis circadian of Botany, 99, 529–536. clock component PRR7. The Plant Journal, 76, 101–114. https​://doi. Savolainen, O., Pyhäjärvi, T., & Knürr, T. (2007). Gene flow and local org/10.1111/tpj.12276​ adaptation in trees. Annual Review of Ecology, Evolution, and Liu, X., & Fu, Y.‐X. (2015). Exploring population size changes using Systematics, 38, 595–619. https​://doi.org/10.1146/annur​ev.ecols​ SNP frequency spectra. Nature Genetics, 47, 555–559. https​://doi. ys.38.091206.095646 org/10.1038/ng.3254 Schwaegerle, K. E., Mcintyre, H., & Swingley, C. (2000). Quantitative Lotterhos, K. E., & Whitlock, M. C. (2014). Evaluation of demographic his‐ genetics and the persistence of environmental effects in clon‐ tory and neutral parameterization on the performance of FST outlier ally propagated organisms. Evolution, 54, 452–461. https​://doi. tests. Molecular Ecology, 23, 2178–2192. org/10.1111/j.0014-3820.2000.tb000​48.x Martin, M. (2011). Cutadapt removes adapter sequences from high‐ Tajima, F. (1989). Statistical method for testing the neutral mutation hy‐ throughput sequencing reads. EMBnet.journal, 17(1), 10–12. https​:// pothesis by DNA polymorphism. Genetics, 123, 585–595. doi.org/10.14806/​ej.17.1.200 Thomas, P. A., & Polwart, A. (2003). Biological flora of the British Isles. Mayol, M., Riba, M., Cavers, S., Grivet, D., Vincenot, L., Cattonaro, F., … Journal of Ecology, 91, 489–524. González‐Martínez, S. C. (2019). Data from: A multiscale approach Tiffin, P., & Ross‐Ibarra, J. (2014). Advances and limits of using population to detect selection in non‐model tree species: widespread adap‐ genetics to understand local adaptation. Trends in Ecology & Evolution, tation despite population decline in Taxus baccata L. Dryad Digital 29, 673–680. https​://doi.org/10.1016/j.tree.2014.10.004 Repository, https​://doi.org/10.5061/dryad.78j501j Tintle, N. L., Borchers, B., Brown, M., & Bekmetjev, A. (2009). Comparing Mayol, M., Riba, M., González‐Martínez, S. C., Bagnoli, F., de Beaulieu, gene set analysis methods on single‐nucleotide polymorphism data J. L., Berganzo, E., … Vendramin, G. G. (2015). Adapting through from Genetic Analysis Workshop 16. BMC Proceedings, 3, S96. https​ glacial cycles: Insights from a long‐lived tree (Taxus baccata). New ://doi.org/10.1186/1753-6561-3-S7-S96 Phytologist, 208, 973–986. Turner, C. (2000). The Eemian interglacial in the North European plain Mwimba, M., Karapetyan, S., Liu, L., Marqués, J., McGinnis, E. and adjacent areas. Netherlands Journal of Geosciences, 79, 217–231. N., Buchler, N. E., & Dong, X. (2018). Daily humidity oscilla‐ https​://doi.org/10.1017/S0016​77460​0 023660 tion regulates the circadian clock to influence plant physiol‐ Villemereuil, P., Frichot, E., Bazin, E., Francois, O., & Gaggiotti, O. E. ogy. Nature Communications, 9, 4290. https​://doi.org/10.1038/ (2014). Genome scan methods against more complex models: When s41467-018-06692-2 and how much should we trust them? Molecular Ecology, 23, 2006– Nakamichi, N., Kita, M., Niinuma, K., Ito, S., Yamashino, T., Mizoguchi, T., 2019. https​://doi.org/10.1111/mec.12705​ & Mizuno, T. (2007). Arabidopsis clock‐associated Pseudo‐Response Villemereuil, P., & Gaggiotti, O. E. (2015). A new FST‐based method to Regulators PRR9, PRR7 and PRR5 coordinately and positively reg‐ uncover local adaptation using environmental variables. Methods in ulate flowering time through the canonical CONSTANS‐dependent Ecology and Evolution, 6, 1248–1258. photoperiodic pathway. Plant and Cell Physiology, 48, 822–832. https​ Villemereuil, P., Gaggiotti, O. E., Mouterde, M., & Till‐Bottraud, I. (2016). ://doi.org/10.1093/pcp/pcm056 Common garden experiments in the genomic era: New perspectives Nazareno, A. G., Bemmels, J. B., Dick, C. W., & Lohmann, L. G. (2017). and opportunities. Heredity, 116(249), 254. https​://doi.org/10.1038/ Minimum sample sizes for population genomics: An empirical study hdy.2015.93 160 | MAYOL et al. Whitlock, M. C. (2000). Fixation of new alleles and the extinc‐ Yeaman, S., Hodgins, K. A., Lotterhos, K. E., Suren, H., Nadeau, S., tion of small populations: Drift load, beneficial alleles, and Degner, J. C., … Aitken, S. N. (2016). Convergent local adaptation to sexual selection. Evolution, 54, 1855–1861. https​://doi. climate in distantly related conifers. Science, 353, 1431–1433. https​ org/10.1111/j.0014-3820.2000.tb012​32.x ://doi.org/10.1126/scien​ce.aaf7812 Willi, Y., Van Buskirk, J., Schmid, B., & Fischer, M. (2007). Genetic iso‐ lation of fragmented populations is exacerbated by drift and se‐ lection. Journal of Evolutionary Biology, 20, 534–542. https​://doi. org/10.1111/j.1420-9101.2006.01263.x S U P P O R T I N G I N FO R M AT I O N Willing, E. M., Dreyer, C., & Van Oosterhout, C. (2012). Estimates of ge‐ netic differentiation measured by FST do not necessarily require large Additional supporting information may be found online in the sample sizes when using many SNP markers. PLoS ONE, 7, e42649. Supporting Information section at the end of the article.   https​://doi.org/10.1371/journ​al.pone.0042649 Wright, S. (1931). Evolution in Mendelian populations. Genetics, 16, 97–159. How to cite this article: Mayol M, Riba M, Cavers S, et al. A Wright, S. I., & Andolfatto, P. (2008). The impact of natural selection on the genome: Emerging patterns in Drosophila and Arabidopsis. Annual multiscale approach to detect selection in nonmodel tree Review of Ecology, Evolution, and Systematics, 39, 193–213. species: Widespread adaptation despite population decline in Yang, J., Benyamin, B., McEvoy, B. P., Gordon, S., Henders, A. K., Nyholt, Taxus baccata L. Evol Appl. 2020;13:143–160. https​://doi. D. R., … Visscher, P. M. (2010). Common SNPs explain a large pro‐ org/10.1111/eva.12838​ portion of the heritability for human height. Nature Genetics, 42, 565–569. https​://doi.org/10.1038/ng.608