Aim: Climate change poses significant challenges for tree species, which are slow to adapt and migrate. Insight into genetic and phenotypic variation under current landscape conditions can be used to gauge persistence potential to future conditions and determine conservation priorities, but landscape effects have been minimally tested in trees. Here, we use Pinus contorta, one of the most widely-distributed conifers in North America, to evaluate the influence of landscape heterogeneity on genetic structure as well as the magnitude of local adaptation versus phenotypic plasticity in a widespread tree species.

Location: Western North America.

Methods: We paired landscape genetics with fully reciprocal in situ common gardens to evaluate landscape influence on neutral and adaptive variation across all subspecies of P. contorta.

Results: Landscape barriers alone play a minor role in limiting gene flow, creating marginal geographically-based structure. Local climate determines population performance, with survival highest at home but growth greatest in mild climates (e.g., warm, wet). Survival of two of the three populations tested was consistent with patterns of local adaptation documented for P. contorta, while growth was indicative of plasticity for populations grown under novel conditions and suggesting that some populations are not currently occupying their climatic optimum.

Main Conclusions: Our findings provide insight into the role of the landscape in shaping population genetic structure in a widespread tree species as well as the potential response of local populations to novel conditions, knowledge critical to understanding how widely distributed species may respond to climate change. Geographically-based genetic structure and reduced survival under water-limited conditions may make some populations of widespread tree species more vulnerable to local maladaptation and extirpation. However, genetically diverse and phenotypically plastic populations of widespread trees, such as many of the P. contorta populations sampled and tested here, likely possess high persistence potential.

Methods are extensive, so please refer to manuscript for details.

Twenty sampling locations were randomly selected from occurrences in regions one through four. In region four, we included two sampling locations representing proposed variety yukonensis (Strong, 2010), and we avoided sampling across much of central and southern British Columbia where commercial plantations of latifolia are common. Some region five samples of latifolia were provided by the FIA program (n = 13), with additional sampling locations selected to fill in gaps not sampled by FIA (n = 3). In region six, the two known populations were sampled. Ultimately, fifty-one locations were sampled. At each sampling location, one gram of current-year needles was collected from ten individuals (>50m apart) and preserved using silica gel desiccant. Total genomic DNA was extracted using DNeasy plant kits (Qiagen, Valencia, CA) at the U.S. Department of Agriculture National Forest Genetics Laboratory (Placerville, CA). Of 15 highly polymorphic SSR markers initially tested (Lesser et al. 2012), nine amplified across all samples (Appendix Table A1). Loci were amplified in multiplex under identical conditions, with locus-specific primers 5’-tailed with universal primer sequences (as described by Missiaggia & Grattapaglia 2006, see Appendix 1 for details). PCR products were separated on a 3730xl Genetic Analyzer (Life Technologies, Carlsbad, CA), and peak sizes were determined using GeneMarker v2.2 (SoftGenetics LLC, State College, PA). Samples were scored three times to verify peaks and resolve conflicts.

After screening and adjusting for null alleles, genotyping errors, and deviations from Hardy-Weinberg Equilibrium, we calculated pairwise F_ST (i.e., the inbreeding coefficient or proportion of genetic variance contained within a subpopulation relative to total genetic variance) and the following parameters, averaged across loci, for each sampling location using GenAlEx (Peakall & Smouse, 2012): percent polymorphic loci (PPL), allelic richness (N_A), number of effective alleles (N_E), expected heterozygosity (H_E), unbiased expected heterozygosity (uH_E), and inbreeding levels (F_IS). We used the ‘pegas’ package in R (Paradis, 2010; R Core Team, 2019) to quantify population differentiation within and among sampling locations and subspecies using a hierarchical analysis of molecular variance (AMOVA).

Pairwise genetic distances among sampling locations were calculated using conditional genetic distance (cGD), where genetic distances are based on genetic covariance and estimated from graph distances as the shortest path connecting pairs via population graph topology (Dyer et al., 2010). Pairwise cGD is more sensitive than traditional metrics (e.g., F_ST), accounting for both direct and indirect connectivity (Dyer et al., 2010). We estimated cGD using the ‘GStudio’ package in R (Dyer, 2016).

We tested for range-wide genetic connectivity by comparing pairwise cGD to pairwise spatial and environmental distances, testing hypotheses of isolation by distance (IBD), barrier (IBB), resistance (IBR), and environment (IBE). For tests of IBD, we calculated pairwise Euclidean geographic distance (km) using Vincenty ellipsoid distance in the ‘geosphere’ package in R (Hijmans et al., 2019). For IBB, we created a binary matrix representing hypothesized barriers to gene flow: Central Valley of California separating coastal and mountain populations, Juneau Icefield and Coast Mountain Range separating coastal Alaska and interior Yukon and British Columbia populations, and Great Basin-Intermountain West separating Sierra Nevada and Rocky Mountain populations. Tests of IBR were performed using a resistance map derived from habitat suitability modeling, representing landscape resistance to movement among populations. IBE was evaluated using among-population climate dissimilarities irrespective of spatial connectivity, calculated as pairwise Euclidean distances based on the first three principal components from an analysis of seven bioclimatic variables (‘prcomp’ function in R). We used multiple approaches to evaluate which hypotheses (IBD, IBB, IBR, IBE) best describe observed patterns of genetic distance. First, we used Mantel and partial Mantel tests in the R package ‘vegan’ (Oksanen et al., 2018) under a reciprocal causal modeling framework (Cushman et al., 2013) to evaluate relative support as the difference between reciprocal partial Mantel tests for each hypothesis. Because Mantel and partial Mantel tests are criticized for their tendencies toward inflated type I error rates (Guillot & Rousset, 2013), we also implemented multiple matrix regression with randomization (MMRR, Wang, 2013) in the R package ‘ecodist’ (Goslee & Urban, 2007) to test for consistency of results, comparing all possible combinations of hypotheses to identify the models with the greatest support.

Note that Population 23 (contorta9 in manuscript) was only incldued in population genetics summary in GenAlEx. Samples from contorta 9 (Jumbo Bog, Alaska) were difficult to score at numerous loci; we included this sampling location in population frequency statistics but excluded it from all population and landscape genetics analysis.

All climate data was extracted from TerraClim.