Conservation prioritization based on past cascading climatic effects on genetic diversity and population size dynamics: Insights from a temperate tree species

The biogeographic history affects adaptive evolution by altering the pattern of genetic variation and provides the background upon which other evolutionary processes are operating. Assessment of the consequence of such eco‐evolutionary interactions is crucial to population's genetic diversity maintenance. Yet, considering how biogeographic and multiple evolutionary processes interact in complex environments to alter genetic diversity in conservation planning remains unresolved.


| INTRODUC TI ON
Anthropogenic climate change is causing a wide-spread maladaptation in natural ecosystems and poses a potent threat to global biodiversity (IPBES, 2019;Shaw & Etterson, 2012;Urban, 2015).
Forecasting the loss of suitable habitats caused by climate change helps inform conservation interventions and provide guidelines and tools for assisted adaptation and migration (McLachlan et al., 2007;Thomas et al., 2004;Thuiller et al., 2005). Models currently used to predict species responses to climate change often utilize a correlative niche modelling approach (Thuiller et al., 2008), which assumes static environmental tolerance over time such that organisms will track their current ecological niche in the future (Araújo & Peterson, 2012). However, this type of models could spuriously forecast less future suitable habitats because species' current distributions may not adequately represent their thermal tolerance, and future climate conditions may not be analogous to current conditions (Williams & Jackson, 2007). Decades of empirical research also contradict this assumption, as evolutionary processes are also able to counteract the impact of climate change in some species (Franks et al., 2007;Karell et al., 2011), while in others, the strong selective pressures imposed by climate change might outpace a species' capability of adapting to a changing environment (Etterson & Shaw, 2001;Jezkova & Wiens, 2016;Quintero & Wiens, 2013). In fact, the feedback between adaptive and demographic processes has large consequences for a species' response to climate change (Hoffmann & Sgrò, 2011;Urban, 2015). Hence, examination of demographic, ecological and evolutionary mechanisms influencing population divergence improves our prediction of the impacts of climate change on ecologically divergent populations.
Historically, past glaciation, most notably the Pleistocene ice ages, has profoundly shaped the genetic architecture of the flora and fauna of the Pacific Northwest and created discontinuities in the geographic distribution of many plant species (Milne, 2006;Soltis et al., 1997). Such spatial patterns of glacial persistence and post-glacial colonization led to dramatic demographic shifts of populations or species, and theoretically, large populations recover to equilibrium diversity levels more slowly (Kimura, 1983 p. 204;Nei & Graur, 1984 p. 82-89). Moreover, past demography determines the strength of genetic drift [i.e. random changes in allele frequency over time; (Caballero, 1994)] and has a profound influence on efficiency of selection (Wright, 1931) as well as genetic diversity (Hewitt, 2000) such as possible genetic erosion after a strong bottleneck. In addition, fluctuations in population size due to historical dynamics of population abundance are assessed widely as an important parameter for deriving neutral expectations (Lotterhos & Whitlock, 2014). Hence, reconstructing the demographic history of populations or species is paramount to understanding environmentdriven population divergence and providing useful information for conservation planning.
Long-lived organisms, such as trees, are particularly vulnerable to maladaptation due to climate change, because the long life span limits their adaptive capacity (Aitken et al., 2008;Kuparinen et al., 2010).
Cottonwoods (Populus spp.) are a genus of deciduous flowering trees, in which black cottonwood (Populus trichocarpa) represents an important model species for the long-lived perennial plants (Box 1). Genetic analysis offers evidence that post-glacial recolonization from multiple refugia (e.g. one or more northern refugium) is a plausible scenario for the demographic history of P. trichocarpa Zhou et al., 2014). In addition, a glacial refugium of Populus trees has been located to the panhandle of Alaska or Beringia, an unglaciated region encompassing the former Bering land bridge and the land between the Lena and Mackenzie rivers (Breen et al., 2012;Brubaker et al., 2005;Levsen et al., 2012). At present, P. trichocarpa inhabits riparian areas in western North America from Baja California to Alaska and grows preferably west of the Cascade Mountains in northern Oregon, Washington, and southern British Columbia (DeBell, 1990).
This riparian species is highly sensitive to changes in water levels (Rood et al., 2003;Rood & Mahoney, 1990) which climate change can affect by altering snowpack melt and rainfall.

BOX 1 Use of Populus trichocarpa as a model system for genetic and evolutionary studies
Populus trichocarpa has porous species barriers, moderate genome size (c. 480 Mb) and abundant genomic resources, rendering it an ideal system for molecular evolutionary studies (Cronk, 2005;Jansson & Douglas, 2007;Tuskan et al., 2006). The genomic basis of adaptive variation across a wide geographic range of P. trichocarpa individuals has been elucidated (Evans et al., 2014), and the adaptive and neutral evolutionary forces shaping the pattern of P. trichocarpa population structure comprise introgression and hybridization with the congener P. balsamifera (Suarez-Gonzalez et al., 2018), restricted gene flow (Geraldes et al., 2014;Xie et al., 2009), natural selection (Geraldes et al., 2014) and genetic drift (population history) . Extensive gene flow is facilitated by long-distance dispersal through pollen, seeds and vegetative propagules (Galloway & Worrall, 1979) with effective distances of 2-7.6 km Slavov et al., 2009). In addition, phenotypic plasticity has been shown to be another important factor affecting trait variability and adaptability (Liu & El-Kassaby, 2019). At the genomic level, divergent selection is evident in promoters, 5'-UTR and intronic sites, whereas intergenic regions evolve under relaxed selection (Zhou et al., 2014). Genetic divergence genes of introgression-resistant regions are found to be under positive selection (Liu & El-Kassaby, 2021). A single-genebased mechanism determines the sexes of this dioecious species, whereby impacting species evolution (Müller et al., 2020). Latitudinal allele frequency clines are also common across the genome, hinting at large-scale patterns of adaptation in the species .
Previous studies suggest that effective population size (N e ) P. trichocarpa has ranged between 4000 and 8000 individuals (Levsen et al., 2012;Slavov et al., 2012;Zhou et al., 2014); however, this estimate is remarkably low with respect to the vast, current census population size of P. trichocarpa stands and three to four times lower than those that would be obtained based on putative neutral nucleotide diversity (Gilchrist et al., 2006;Tuskan et al., 2006). Moreover, previous work reported the highest N e in the northern population (Zhou et al., 2014), which is counterintuitive given that the presentday core growth region is in southern British Columbia and northern Oregon.
To understand the way the P. trichocarpa populations respond to climate change depends on the extent to which the biogeography history has shaped their genetic divergence, I use multifaceted approaches to address the question of how past climate has affected the dynamics of niche evolution, demographic history and genetic diversity. To that end, I specifically: • decipher population structure to identify genetic clusters, followed by characterizing distinctiveness via genomic metrics (e.g. π, F ST and number of private alleles) between clusters ("populations" thereafter) and assume each population with its own unique evolutionary trajectory; • scan for recent positive or long-term balancing selection in the genome of each population and assume that populations with more adaptive genetic loci indicate their great evolutionary potential because selective sweeps could lead to shifts in the selection landscape and thus local adaptation; • estimate the direction and magnitude of gene flow between populations, whereby inferring the demographic history for single populations and the entire species; • perform ecological niche modelling to forecast habitat suitability under climatic scenarios and test for niche evolution (divergence or conservatism); and • predict changes in population genetic structure in response to climate change.
By assembling the aforementioned information (Figure 1), I provide the idea that we could facilitate prioritization of target and source populations for proactive management interventions through assisted gene flow to prevent populations from genetic loss.

| Plant material and georeferenced data
The Populus trichocarpa provenances used in this study were collected by British Columbia Ministry of Forests, Lands and Natural Resource Operations (Xie et al., 2009), and planted in a common garden at Totem Field of the University of British Columbia. A total of 409 individuals from 139 provenances were selected. These individuals were from 29 drainages (i.e. topographic units separated by watershed barriers) extending 14° in latitude (45. 6°N-59.6°N) within its geographic distribution range (31°N-62.5°N) (Appendix S1: Table A1). Georeferenced occurrence records for P. trichocarpa were retrieved from the literature and herbaria accessed through the Global Biodiversity Information Facility data portal at https:// www.gbif.org/. A total of 886 distinct locales were gleaned and used to represent the entire range of P. trichocarpa (Appendix S1: Table   A2). Samples were grouped into three present-day populations including British Columbia (BC)-North, BC-South and Oregon.

| Genotype data
After the raw reads were filtered into clean biallelic SNPs using a set of filtering criteria (Method S1 and Figure S1), a total of 1.6 million SNPs (Dataset I) were extracted. Population parameters (Nei's nucleotide diversity (π), Tajima's D, linkage disequilibrium (LD) statistics, Wall's B and Z nS and Hudson's F ST ) were calculated in a 500-bp sliding window with a step size of 500-bp on PopGenome v.2.2.4 (Pfeifer et al., 2014). The diversity values were normalized by the number of total nucleotides (both variable and invariable sites) in each window. Based on the SNPs previously identified through associations with phenotypes and/or habitat climates (Table S1 and Figure S2), I aggregated 681 SNPs, termed genome-wide associated (GWA) SNPs (Dataset II), for subsequent use (Appendix S1: Table A3).
To decrease false positives, all retained SNPs passed the filtering of being identifiable in two versions of the P. trichocarpa genome (v.2.0 and 3.0; downloaded from Ensembl Plants). Likewise, the population parameters were calculated for each population in a 1-bp sliding window with a step size of 1-bp using PopGenome v.2.2.4 (Pfeifer et al., 2014). Hardy-Weinberg equilibrium (HWE) and LD were calculated using VCFtools v.4.2 (Danecek et al., 2011). The false discovery rate (Benjamini & Hochberg, 1995) was applied to correct P-values from HWE and LD multiple exact tests with α = 0.05 using the q value package v.2.14.0 (Bass et al., 2015) in R. As the combined effects of the Earth's tilt and rotation generate annual sine wave variation in day length with cascading effects on temperature, rainfall and resource availability, the average day length | 2715 LIU (lenDay; h) in June (summer solstice included) was calculated for each locale as a proxy for the photoperiodic regime using the geosphere package v.1.5.10 in R. I further extracted 19 bioclim variables from WorldClim v.2.0 (only available for the present-day conditions; last access on April 2020) and v.1.4 (past and future conditions) (ref. Hijmans et al., 2005; http://world clim.org) at 2.5 arc-minutes resolution. These data allowed for modelling a species' full distribution range only via latitude and longitude extraction from a Spatial Polygons Data Frame or shapefile. A total of 15,602 spatial points within its natural range were extracted from a digitized map in ref. (Little, 1971), accessible from Data Basin at https://datab asin.org/.

| Climate data
As the existence of strongly correlated variables may overfit the model, I reduced collinearity (|r| > 0.

| Population genetic assignment
To understand population genetic structure, I resorted to multiple analyses including SNP-, frequency-, distance-and likelihood (Bayesian)-based approaches. Multiple genetic differentiation metrics were calculated for populations using PEAS v1.0 (Xu et al., 2010) (Method S2). Compared to 34K genome-wide markers (Geraldes et al., 2014) previously used for genetic structure analysis, I used neutral GWA SNPs (N = 577 based on Dataset II; Method S3, Figures   S4 and S5 for their identification details). Principle component F I G U R E 1 Conceptual representation of multiple genomic analyses to identify donor/recipient populations for developing conservation strategies. Illustrations are to elucidate different approaches employed to compare populations (e.g. Popltn 1° and 2°) based on (1) prior knowledge on eco-evolutionary interactions including genomic metrics/diversity (coloured dots), signals of selection across the genome (Manhattan plot), and inferred demographic histories (changes in population sizes through time); and (2) predictions of change under future climatic scenarios using species niche suitability (shrinking or expanding geographic distribution ranges) and genetic makeup turnover (genetic structure differentiated by colour) great genetic turnover Population genetic assignment High → low diversity π; Great → low evolutionary potential (e.g., by # of adaptive loci identified as an indicator); Large → small Ne (e.g., conserve populations with recent bottleneck); High → low niche suitability (e.g., expanded → shrunk population); Small → high changes in population genetic structure.  (Lȇ et al., 2008) in R.

Assisted migration or gene flow
As the PCA does not entail the consideration of population structure, HWE and LD, I further implemented a Bayesian approach in STRUCTURE v.2.3 to infer population structure (Pritchard et al., 2000). This program iteratively ascribes individuals to K panmictic groups (populations) by minimizing deviations from HWE and LD within groups. I assumed that individuals followed the admixture model in structure, and allele frequencies were correlated among groups (Falush et al., 2003). For the correlated allele frequency model, the parameter λ was estimated at K = 1, as recommended by the software documentation for SNP data with rare alleles.
According to a previous empirical estimate of the number of genetic groups using a Maximum Likelihood model (Geraldes et al., 2014), I varied K from 1 to 15, and for each run, I conducted five replicates with a burn-in period of 10,000 and MCMC iterations of 100,000.
I applied posterior log probability (Pritchard et al., 2000), [Pr(X|K)], for each K and ad hoc statistic (Evanno et al., 2005), ΔK, to calculate the uppermost hierarchical level of structure (K). The results from five independent runs were summarized, and ln Pr(X|K) and ΔK were, respectively, plotted against K using the program STRUCTURE HARVESTER (Earl & vonHoldt, 2012). I used the computer program CLUMPP v1.1.2 to analyse the results from replicate analyses for optimal alignments of replicate clusters (Jakobsson & Rosenberg, 2007). The output from CLUMPP was graphically displayed by the cluster visualization program DISTRUCT v1.1 (Rosenberg, 2004).
Due to significant patterns of IBD ( Figure S7) and geographic barriers ( Figure 2a) between populations, I also used conStruct v1.0.3 (Bradburd et al., 2018) to examine fine-scale population structure in a spatially aware context. Conceptually, conStruct is akin to STRUCTURE, except that geographic distances between sampling locales are included in the model and used for population structure inference. Due to a small number of SNPs used (N = 577), it was not possible to use a cross-validation to find the best-fit number of genetic clusters through conStruct. Instead, a set of spatial models was run by using different cluster numbers (K = 2 to 5 layers), each model with four independent MCMCs and 50,000 generations per MCMC.
Convergence of independent analyses was assessed based on trace plots generated by conStruct.
Furthermore, I tested whether population structure is related to IBD or isolation-by-environment (IBE). Specifically, I tested for significant correlations between matrices of Euclidean genetic, geo- were used to test differentiation strength between groups.

| Genome-wide scan for selective sweeps
Scans for selective sweeps were carried out first using integrated haplotype scores (iHS) implemented in the rehh package v.2.0.2 (Gautier et al., 2017) in R. The rationale for this method is that increase in the frequency of a beneficial allele occurs too quickly for recombination to homogenize its genetic background. Hence, a positively selected allele is embedded in a less common haplotype that conducive to the maintenance of standing genetic variation, by estimating normalized β statistics implemented in BetaScan (Siewert & Voight, 2017). This measure identifies the clusters of variants with an excess number of intermediate frequency polymorphisms (Siewert & Voight, 2017). I used a window size of 1 kb to calculate the unfolded version of β (allelic state as described above) for each core SNPs in each of the three populations. SNPs with extreme β scores in the top 0.05% were deemed significant.

| Historical and contemporary gene flow
To investigate gene flow scenarios between population groups, I employed MIGRATE-N v3.6.11 (Beerli, 2006;Beerli & Felsenstein, 2001;Beerli & Palczewski, 2010) and BayesAss v3.0.4 (Beerli & Felsenstein, 2001), which estimate historical and contemporary migration rates, respectively. I used a Bayesian method to estimate historical gene flow parameters: mutation-scaled effective population where N e is effective population size and μ is the mutation rate per generation) and mutation-scaled migration rates, The program calculates unidirectional estimates of m for each group pair (Fraser et al., 2007), and this method is based on the model in which individuals are exchanged between groups over generations (Goossens et al., 2005). I run the program with the common-line flag "BA3 -v -i10000000 -b1000000 -n1000" on BayesAss (Beerli & Felsenstein, 2001). As the data exceeded the program's single load size, I divided the data into two equal sets (i.e. 341 or 340 loci per load), and the averaged output was used for the final results.
To test scenarios of gene flow among provenances of the three populations, I carried out a multivariate provenance (subpopulation) graphing approach via PopGraph v.1.4 (Dyer & Nason, 2004) using gstudio (Dyer, 2009

| Demographic history inference
The single population demographic histories were estimated using 200 bootstrap iterations with stairway plot package v.2.0 beta2 (Liu & Fu, 2015) based on observed unfolded SFS from each population.
Considering that the exclusion of singletons would skew the entire frequency spectrum potentially impacting population genetic processes inference, I included singletons when estimating the SFS. I used default 2/3 of the data for training, and four numbers of sequences (nseq), including (nseq)-2)/4, (nseq-2)/2, 3*(nseq-2)/4 and nseq-2 as the number of random breakpoints. According to the stairway plot pipeline, the best number of random breakpoints was chosen based on the training data. To calibrate the results, I used a generation time of 15 years (Ingvarsson, 2008) and two different mutation rates, μ = 1.75 × 10 −8 and 2.5 × 10 −9 (Method S8).
To elucidate the population split history, I used the joint 2D-SFS (Method S9) based on the composite likelihood (CL) demographic modelling approach described by Excoffier et al. (2013) and implemented in fastsimcoal v.2.6. I assumed two models with alternative scenarios of lineage divergence or formation: Model one assumes migrants from BC-South (source) to BC-North (sink) due to admixture or introgression events, whereas Model two assumes the same events with a reverse direction of migration. In both models, population sizes were allowed to vary in historical events, and thus the bottleneck may occur given new population size inferior to the previous one in simulations. The historical gene flow was set according to the result shown in Figure S13. To enable the estimation of other parameters (Excoffier et al., 2013), the parameter for ancestral population size (N e ) was anchored to be 3048, estimated by the formula, N e = N ANC = θ π × (4 × μ × L) −1 , where the genetic diversity (θ π ) was set as 0.0032 (Levsen et al., 2012), the neutral mutation rate (μ) was set as 1.75 × 10 −8 (Method S8), and the generation time (L) was assumed to be 15 years (Ingvarsson, 2008). The recombination rate per generation (c) was set as 9.37 × 10 −8 given c = ρ × (4N e ) −1 , where ρ is the scaled recombination rate (Zhou et al., 2014). Note that these two N e estimation approaches may not be comparable, as SFS data and short-range LD are very sensitive to long-term N e , whereas long-range LD tends to be influenced by very recent changes in population size (e.g. bottleneck and recent expansion) (Rogers, 2014;Tenesa et al., 2007).

| Ecological niche modelling and niche divergence
To test how climate has driven species adaptation, I performed climatic suitability prediction using ecological niche modelling ( Predicted distribution probability (on a logistic scale) in each grid cell is an index of suitability or suitability scores, ranging from 0 to 1 (unsuitable to highly suitable [blue to red portrayed on the graph]).
I produced a current distribution using default settings except only allowing hinge features, a 20-fold cross-validation, and a regularization multiplier of 2.5 to fit more general models (Elith et al., 2010). MCMC runtime was set to 100,000 sweeps and the burn-in period to 10,000 sweeps, and the degree of the spatial trend surface was set to 1. These settings were launched using models with admixture.
To define the number of clusters, K, I run the simulations five times | 2719 LIU for each K ranging between 2 and 15. A subset of runs that minimize the deviance information criterion (DIC) scores (Spiegelhalter et al., 2002), namely, the lowest DIC values, were selected. Finally, I forecasted the genetic clusters for two future scenarios (RCP 4.5 and 8.5) under two GCMs: CCSM4 and MIROC. To average multiple runs with the same number of clusters, I used the computer program CLUMPP to analyse the results from replicate analyses for optimal alignments of replicate clusters (Jakobsson & Rosenberg, 2007). The output from CLUMPP was graphically portrayed using R. The grid from an ascii raster file created in P. trichocarpa distribution modelling under the current scenario was used as background to display estimated and predicted coefficients of this analysis on the map.

| Population structure revealed by associative genomic loci in a spatially aware context
According to a previous genetic clustering based on genome-wide

| Contribution of IBD and IBE to population genetic divergence
Genetic differentiation between BC-South and Oregon was highest (Wright's F ST = 0.16) (Table S3). Molecular variance by AMOVA (Table   S4) revealed that most of genetic variation was distributed within individuals (86.04%), and individuals were genetically more divergent among populations (F IT = 0.14) than within population (F IS = 0.06).
The spatial autocorrelation analysis showed significantly positive genetic correlations for populations within 200 km ( Figure S7). There was a moderately high correlation between climate and geography (r = .627, p < .0001) ( Figure S8a). The correlation between genetic and climatic distances increased (r = .49, p < .0001) when the three populations were pooled for analysis ( Figure S8b). Note that drainages #29 and #30 of Oregon formed two distinct clouds of points, and the correlation was weak based on either drainage (average r = .23, p < .0001) ( Figure S8b). Nonetheless, some climatic variables (e.g. mean annual precipitation, MAP) were weakly-to-moderately correlated with niche climate and geography (r = .13-.49, p < .0001) ( Figure S9). In addition, average day length (photoperiod) was highly correlated with geography as expected, showing a latitudinal pattern (spearman's ρ = 0.955, p < .0001; Figure S2), but the correlation between day length and niche climate was much weaker (ρ = 0.629, p < .0001; Figure S2). There were intermediate-to-weak correlations of MAP and mean annual temperature (MAT) with geography (ρ = 0.594-0.132, all p < .0001; Figure S9).

| Ecological niche model and niche evolution reveal niche suitability and divergence over time
Six of the original 19 bioclimatic variables were selected for ecological niche modelling (ENM) based on collinearity among the variables (Table S2 and Figure S3). The ENM accuracy was relatively high southern Alaska), relative to the present ( Figure S10). I additionally F I G U R E 3 Genome-wide scan for genetic diversity, LD, positive and balancing selection in each population. (a) A 500-bp window size with a 500-bp step size was used to define the slide window for scanning across the 19 chromosomes using the filtered sequentially concatenated 1.6 M SNPs. The π values were weighted by using the total genome length (i.e. both variable and invariable sites). Wilcoxon signed-rank tests were performed between each population pairs. All p-values were adjusted using the sequential Bonferroni method (***:

| Direction of historical and contemporary gene flow between populations
Historical gene flow was low and unidirectional from BC-North to -South, but there was a strong and bidirectional gene flow between BC-South and Oregon ( Figure S13). Contemporary gene flow within populations was higher than between populations ( Figure   S13), and this was also supported by high genetic connectivity among provenances within populations ( Figure S14a). However, high F I G U R E 4 Ecological niche model, niche evolution and demographic inference. (a) Changes in niche suitability for 524 P. trichocarpa locales between present-day  and Last Glacial Maximum (LGM) based on two models, CCSM4 and MIROC. Hindcasted suitability maps can be viewed in Figure S10. An asterisk * for significance at p < .05 based on Wilcoxon test. Non-significant pairs are not displayed. (b) Assessment of niche divergence and conservatism based on pairwise niche divergence tests (Schoener's D). Histograms for null distributions indicate the background levels of niche divergence, and dashed lines indicate the observed value of Schoener's D (niche overlap) for each pair compared. An overlap value larger than the null distributions indicates niche conservatism, whereas a value significantly smaller than the distributions indicates niche divergence. All p-values from background tests are smaller than 0.01. (c) Estimated single population demographic histories. The bold line represents the median estimate from 200 bootstrap replicates with 95% CIs shown via light grey shading. MIS4 denotes Marine Isotope Stage 4. Mutation rate was set as 1.75 × 10 −8 per site per generation, and generation time was assumed to be 15 years in the model. Other demographic trajectories modelled with 7× lower mutation rate can be viewed in Figure S15. (d) Schematic illustration of the best-fitting demographic model for the P. trichocarpa colonization history. The width of the boxes represents the relative effective population size and arrows represent directional gene flow between populations. Maximum-likelihood parameter estimates for the model are given in the graph, including split times (Ka), population sizes and migration rates (95% CIs given in Table S6). The model assumes that BC-North and -South split from a common ancestral and Oregon then separated from BC-South. The three populations underwent recent genetic bottlenecks. Fit of models was assessed and given in Figure S16 and Table S6 BC-North Schoener's D Frequency between-provenance connectivity was not homogenous, and some provenances were dominant gene flow donors or recipients based on low "closeness" and high "degree," "betweenness" and "eigenvector centrality" in PopGraph (Table S5). Moreover, the betweenprovenance correlations were overall moderate within BC-North, -South and Oregon (average Spearman's r = .63, .56, and .35 with average p = .002, .001, and .456, respectively) ( Figure S14b).

| Demographic history shows more recent bottlenecks in BC-North and -South than in Oregon
The stairway plot revealed that bottlenecks in the three popula- The same analysis based on a sevenfold lower mutation rate also revealed bottleneck and gradual N e decline over time in BC-North and -South, respectively, and stable N e in Oregon over the past 60 Ka (comparison with a shorter time period due to a smaller mutation rate used; Figure S15).
Further, I studied the population split times using a simulationbased approach implemented in fastsimcoal2. The model and parameter search ranges were selected based on the results in previous sections (see Methods for setting details; full estimates in Table   S6). Two alternative models with different unidirectional gene flow between BC-North and -South were tested. The best-fitting model favoured unidirectional gene flow from BC-North to -South (Table   S6 and Figure S16). The model also supported the proposed genetic bottlenecks that the species might have undergone (Figure 4d), highlighting more recent bottlenecks in BC-North and -South than in Oregon (Figure 4c). The model also indicates that the north-to-south gene flow into the Oregon was higher than that in the opposite direction (Figure 4d), as shown by one-magnitude higher migration rate for BC-South to Oregon than the reverse direction and 1.4 times higher for BC-North to Oregon than the opposite (Table S6). Under the best-supported model, the current N e of Oregon was largest, in line with the results from the single population demographic models ( Figure 4c; but see Table S7 based on an LD approach using GWA SNPs with a general interpretation in Section 2.7): ~31 million with a 95% CI of 29-52 million from a large ancestral population (~225 million, 95% CI = 98-217 million) due to a bottleneck before separating with BC-South ~51 Ka (95% CI =67-206 Ka) (Figure 4d and Table S6).
BC-North and -South diverged ~610 Ka (95% CI =606-617 Ka) from a large ancestral population (~1.04 million, 95% CI = 0.9-1.05 million) (Figure 4d and Table S6). The simulated SFS from the model showed that both BC-North and -South have an excess of low frequency alleles and a deficit of mid-to-high frequency alleles, consistent with the decline of range-wide empirical SFS ( Figure S17). Such an SFS distribution pattern in the two populations fits with the inference of bottlenecks that they have likely undergone. The Oregon exhibited an exponential decay ( Figure S17a-c) indicating population stationarity, concurring with the outcome of its single population demographic model (Figure 4c).

| Projections of change in niche suitability and genetic structure in response to global warming
Compared to the present, four future distribution models consistently showed no dramatic range shifts over the next 30 years, that is, two generations given 15 years to reproductive maturity in Populus (Δ ≈ ±0.05 or ±5% averaged across the four models; Figure 5a and Figure S11). However, due to the high spatial landscape heterogeneity, findings showed that across populations, niches along the coastal area will likely be less suitable, while the interior has a propensity for increased suitability (Figure 5b). In predicted changes in genetic makeup due to climate change, eight genetic clusters were used based on the deviance information criterion (DIC) scores (Table S8)

| Gene flow and environmental drivers of population divergence
Environmental heterogeneity influences gene flow through geographic barriers that impede dispersal or by creating contrasting selective pressures which elicit local adaptation and select against gene flow between locally adapted forms (Rundle & Nosil, 2005;Sobel et al., 2010). Such influence varies from one system to another, depending on the magnitude and spatial configuration of landscape heterogeneity. In this study, both the best-fitting demographic models (Table S6 and Figure S16 Figure S7). However, there is a climatic similarity between BC-South and Oregon ( Figure S4a), indicating that IBE is not the main cause of their genetic separation. Moreover, IBE was weak within population, as evidenced by significant divergence between adjacent drainages #29 and #30 ( Figure 2d and Figure S8b). This indicates a limited contribution of niche climate adaptation or geographic separation to population structure within P. trichocarpa populations. Furthermore, 11 frost and drought-related climate variables were identified as potential selective drivers ( Figure S5a), in which the primary three included F I G U R E 5 Forecasted changes in niche suitability and genetic clusters due to global warming. (a, b) Changes in niche suitability for 524 P. trichocarpa locales between present-day  and a future period (2050s) based on two models, CCSM4 and MIROC. The future climate predictions considered low and high greenhouse gas scenarios, RCP4.5 and 8.5. Increase (red) and decrease (black) in suitability are respectively graphed in (a) with error bars. The grey-shaded background in (b) is the distribution range of the species. The dashed horizontal lines in (b) mark the boundary for BC-North, Oregon and BC-South, defined as north of 53°N, south of 45°N and between the two transects, respectively. Future suitability maps can be viewed in Figure S10. (c) Spatial distribution of genetic clusters conditioned to climatic variables under low and high greenhouse gas scenarios. Each colour set in the legend represents a genetic cluster. Only grids with suitability above 0.1 projected for the present-day scenario  in the ENM were used for clustering and at each point, the cluster for which the coefficient is maximal was plotted mean annual temperature (MAT), mean warmest month temperature (MWMT) and frost-free period (FFP) ( Figure S5b). However, these climatic stressors were not from single climate clusters ( Figure S4b,c), suggesting that specific climatic variables rather than niche climate similarity (or an aspect of climate) are the true drivers of population genetic divergence.

| Relationship among genetic diversity, selection and demography
Genetic diversity is a key basis for species survival and ecosystem functioning (Toro & Caballero, 2005). Greater genetic diversity prompts higher evolutionary potential and a greater capability of responding to climate change (Frankham et al., 2002). Conserving genetic diversity therefore becomes crucial and high-priority strategies for biodiversity conservation (Frankham et al., 2002).
This study provides evidence that eco-evolutionary processes interact to affect genetic diversity. First, the Tajima's D pattern SFS distribution was graphically used to reveal whether a population is non-bottlenecked at mutation-drift equilibrium or not (Nei et al., 1976). Deviations from an L-shape would indicate a population bottleneck (Luikart et al., 1998). BC-North and -South did not exhibit L-shaped SFS with slight distortions ( Figure S17a Neutral theory predicts that genetic diversity and population size are positively correlated, according to π = 4N e μ where μ is the per site per generation mutation rate. Nonetheless, selection could suppress diversity due to greater impact of linked selection in large populations (Corbett-Detig et al., 2015) possibly driven by higher rates of adaptive substitutions in large populations (Ohta, 1992). A confluence of these viewpoints suggests that selection and demography jointly shape variation in N e and comprise two main factors affecting genetic diversity (Ellegren & Galtier, 2016).
A reduction in N e , a priori, reduces the efficiency of selection, although increased homozygosity could make selection more efficient by exposing recessive mutations. As opposed to the neutral theory, the impact of selection may lead to a negative relationship between π and N e or between estimated heterozygosity (H e ) and N e (Table S7). Indeed, based on demographic models, the Oregon population had the highest N e (Figures 3d and 4c) but the lowest π across the genome ( Figure 3a). Moreover, an early N e drop in the Oregon coincided with the onset of Marine Isotope Stage 4 (~76 Ka; Figure 4c), a period with a rapid increase in global ice volume (Landais et al., 2004). This climatic event occurred 55-64 Ka earlier than the LGM (~21 Ka) or Holocene (~12 Ka), which coincided with a dramatic N e drop in BC-South and -North, respectively ( Figure 4c). Furthermore, demographic history determines the strength of genetic drift, that is, the inverse of N e represents the magnitude of genetic drift and elevated genetic drift results in increased fixation or loss of alleles and therefore low genetic diversity (Wright, 1931). The lowest diversity in the Oregon population might be explained by a possibly extended bottleneck in the Oregon, in contrast with a more recent N e decline in BC-North and -South ( Figure 4c); therefore, BC-North and -South are likely not in equilibrium and thus retain high π.
High genetic diversity prompts more standing genetic variation for selection to act on and thus increases the probability that a population contains genotypes that are favourable under new environmental conditions (Hughes et al., 2008;Reed & Frankham, 2003). Collectively, the key finding of the study supports the impact of demographic dynamics on genetic diversity, and the consequences of this change can be used to infer the vulnerability and persistence of populations to respond under future climatic scenarios. Given that the Oregon population may experience a prolonged bottleneck and climate change would substantially affect its genetic structure, the Oregon may confront a relatively high risk of demographic collapses as climate changes. This evaluation points to conserving the Oregon relative to the other two populations. Broadly, the study underpins that (1) evidence of recent demographic events alone does not suffice to place a population in conservation priority; and (2) genomic diversity estimate, evolutionary potential, future niche suitability and genetic makeup are equally important to knowing whether the population has undergone a genetic diversity loss or even erosion due to recent demographic processes. which this study is based, and to the Editor and two anonymous reviewers for their helpful suggestions.

CO N FLI C T O F I NTE R E S T
The author declares no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Whole-genome shotgun sequencing data and alignment information have been deposited on SRA under the accession, PRJNA276056.