Phenotypic and Genomic Local Adaptation across Latitude and Altitude in Populus trichocarpa

Abstract Local adaptation to climate allows plants to cope with temporally and spatially heterogeneous environments, and parallel phenotypic clines provide a natural experiment to uncover the genomic architecture of adaptation. Though extensive effort has been made to investigate the genomic basis of local adaptation to climate across the latitudinal range of tree species, less is known for altitudinal clines. We used exome capture to genotype 451 Populus trichocarpa genotypes across altitudinal and latitudinal gradients spanning the natural species range, and phenotyped these trees for a variety of adaptive traits in two common gardens. We observed clinal variation in phenotypic traits across the two transects, which indicates climate-driven selection, and coupled gene-based genotype–phenotype and genotype–environment association scans to identify imprints of climatic adaptation on the genome. Although many of the phenotype- and climate-associated genes were unique to one transect, we found evidence of parallelism between latitude and altitude, as well as significant convergence when we compared our outlier genes with those putatively involved in climatic adaptation in two gymnosperm species. These results suggest that not only genomic constraint during adaptation to similar environmental gradients in poplar but also different environmental contexts, spatial scale, and perhaps redundant function among potentially adaptive genes and polymorphisms lead to divergent adaptive architectures.


Introduction
Predicting the evolutionary and ecological outcomes of environmental change requires an understanding of the genomic basis of quantitative, locally adaptive traits in order to assess the magnitude of allele frequency shifts necessary to track novel climatic regimes. If the underlying genotypeenvironment relationships are generalizable among related species that inhabit shared environments, such predictions could be advanced without exhaustive characterization of the genomic architecture of these traits at the species level. Convergent adaptation occurs when distantly related species independently develop similar phenotypes, whereas parallel adaptation refers to populations or closely related lineages evolving repeated patterns of phenotypic divergence in distinct geographical locations (Ralph and Coop 2010;Westram et al. 2014). The formation of phenotypic parallelism relies on standing variation, gene flow, and shared selective regimes (Faria et al. 2014;Westram et al. 2014). When shared standing genetic variation and adaptive loci are transmitted between populations or species, selection may target the same genetic loci to generate similar phenotypes across "replicate" environments (Westram et al. 2014). Factors such as demographic history, effective population size, extent of geographic and genetic separation between populations, and mutation rate could also influence the level of parallelism (Ord and Summers 2015;Kess et al. 2018), and recent theory suggests that nonrandom repeatability of evolutionary trajectories at the genomic level is likely due in part to low redundancy constraining possible paths to adaptation (Yeaman et al. 2018).
Although there are numerous examples of genetic parallelism or convergence for relatively simple traits (e.g., armor plating in sticklebacks [Jones et al. 2012], coat color in mice [Manceau et al. 2010], toxin resistance in snakes [McGlothlin et al. 2016]), less is known about the extent to which shared environmental pressures lead to similar architectures of adaptation for quantitative traits, which encompass the majority of ecologically relevant traits. Range shifts for plant species associated with past climate change have resulted in steep latitudinal and altitudinal gradients that provide natural experiments to test the degree to which parallel environmental clines lead to parallel phenotypic and genomic changes for quantitative, locally adaptive traits (Oubida et al. 2015;Holliday et al. 2016). Functional studies suggest that some of the same genes govern seasonal dormancy, a key adaptive trait, in both angiosperm and gymnosperm trees. For example, downregulation of FLOWRERING LOCUS T leads to bud set in both aspen and spruce trees, the timing of which varies with a tree's latitudinal origins (Bö hlenius et al. 2006;Gyllenstrand et al. 2007). Moreover, population genomic studies have documented parallel or convergent signatures of selection on additional genes involved in local climatic adaptation Yeaman et al. 2016). Despite efforts to understand repeated evolution within and between plant species, questions remain. Is the genetic basis underlying phenotypic adaptation conserved across lineages along parallel ecological clines? Do shared life history characteristics, ecological context, and selection regimes increase the likelihood of parallel genomic architecture of adaptation? Are the same environmental factors driving phenotypic and genetic adaptation?
Populus trichocarpa is a model tree species native to riparian areas of temperate western North America (Gornall and Guy 2007) and is widely distribution from southern Alaska to California. Patterns of variation across the species range for traits, such as bud phenology, biomass and growth, wood biochemistry, and lignin content, suggest that these traits are subject to spatially varying selection (Gornall and Guy 2007;Gailing et al. 2009;McKown, Guy, et al. 2014;Oubida et al. 2015;Wang et al. 2018). Although several studies of phenotypic and genomic adaptation across latitudinal gradients have been reported for Populus species, less attention has been paid to variation along altitudinal clines (Klepsatel et al. 2014;Halbritter et al. 2015). Altitude resembles latitude in forming a gradient of temperature and precipitation but may differ in other climate drivers such as frost, snowfall, evaporative demand, and solar radiation (Thomas 2011). A recent study in P. trichocarpa found a high level of shared loci under divergent selection across latitudinal and altitudinal clines, which is likely a result of parallel environmental selection (Oubida et al. 2015;Holliday et al. 2016). In this study, we aim to further understand parallel adaptation by investigating phenotypic adaptation and climatic selection among P. trichocarpa populations distributed across altitude and latitude. We combine genome-wide association analysis (GWAS) and genotype-environment association analysis to identify key loci underlying adaptive traits or targeted by climate selection in both transects and assess the extent of genetic parallelism across two parallel clines.

Plant Material and Exome Sequencing
Branch cuttings representing 182 provenances spanning 20 of latitude ( fig. 1) were rooted in a greenhouse and four replicates of each genotype were subsequently randomized to one of four blocks and planted in 2012 at two common gardens located at Critz, VA (36.63 N,80.15 W) and Campbell River,British Columbia,Canada (50.06 N,125.32 W). Genomic DNA was extracted from young leaves with the Qiagen DNeasy Plant Mini Kit (Qiagen, Inc, Valencia, CA). Oligonucleotide baits targeting the exome were designed using Agilent SureSelect eArray software (Agilent Technologies, Santa Clara, CA) based on the v2.0 reference P. trichocarpa genome (http://phytozome.jgi.doe.gov; last accessed July 25, 2019). Library preparation and target enrichment procedures were described previously (Zhou and Holliday 2012;Zhou et al. 2014;Holliday et al. 2016;Zhang et al. 2016 Genotyping and Single-Nucleotide Polymorphism Calling Following demultiplexing, trimming of adapter sequence, and removal of low quality reads (Zhou et al. 2014), reads were aligned to the P. trichocarpa genome with the BWA sampe function (Li and Durbin 2009). Indels were realigned using the GATK IndelRealigner function, and duplicate reads marked and removed with Picard MarkDuplicates and GATK DuplicateReadFilter, respectively (DePristo et al. 2011). Single nucleotide polymorphisms (SNPs) were called with the GATK HaplotypeCaller (https://www.broadinstitute.org/ gatk/; last accessed July 25, 2019) (Poplin et al. 2018). Variants were flagged and removed as low quality if they had the following characteristics: low map quality (MQ < 40); high strand bias (FS > 40); differential map quality between reads supporting the reference and alternative alleles (MQRankSum < À12.5); bias between the reference and alternate alleles in the position of alleles within the reads (ReadPosRankSum < À8.0); and low depth of coverage (DP < 5). Missing genotypes on chromosomes 1-19 were imputed using BEAGLE software (Browning BL and Browning SR 2016

Phenotyping of Climate-Related Traits
Height and diameter were measured after all trees set bud in the fall of 2013. The stages of bud set and bud flush were scored weekly until most trees formed terminal bud or until all trees had a fully expanded leaf, respectively (Frewen et al. 2000;Rohde 2002). Bud set timing was calculated as the days elapsed for trees to reach a fully developed bud from January 1, whereas timing of bud flush was expressed as the number of days to the first fully unfolded leaf. Cold hardiness was measured in 2012 by sectioning lateral shoots into $0.5-cm discs and treated samples with temperature of À8, À14, and À20 C (starting temperature: 4 C; decreasing rate 4 C/h) for 1 h before transferring them back to 4 C to thaw. Control samples were maintained at 4 C throughout the experiment. The electrolytic conductivity of the solution was measured and cold injury index calculated after Hannerz et al. (1999), averaged across temperatures for each clone. Finally, after coppicing trees in the Virgina garden in May 2014, we recorded the average height of the regenerated main branches and the number of regenerated branches in March 2015.

Phenotypic Best Linear Unbiased Predictors
We estimated clonal best linear unbiased predictors (BLUPs) for each trait with a linear mixed model using the lmer function in lmerTest package in R (Kuznetsova et al. 2017) as follows: where y ijk is the phenotype of the kth replicate of jth clone in ith block, l is the grand mean, b i is the random effect of ith block, c j is the random effect of jth clone, and e ijk is the random error following Nð0; r 2 e IÞ, where r 2 e is variance and I is the identity matrix. The resulting BLUPs were used as phenotypes in the following GWAS.

Association Analysis of Climate-Related Traits
To control for population structure in our association analyses, we evaluated the fit of several alternative models, namely a simple model (no population structure controlled), models with principal components as covariates (2, 5, and 10 of the leading principal components as covariates), models with familial kinship controlled (K; estimated in TASSEL [Bradbury et al. 2007], and the combined model (2PC þ K), and selected the best structured models using the genomic inflation factor k GC and QQ plots (supplementary table S1 and supplementary fig. S2, Supplementary Material online). With dense SNP coverage and linkage among SNPs, as well as polygenic genetic inheritance underlying the phenotypic traits, the false positive rate is expected to increase (Yang et al. 2011;McKown, Kl ap st e, et al. 2014). As a compromise between false positives and false negatives, we selected the most parsimonious model with a genomic inflation factor (k GC Þ of at least 0.98 but not more than 1.22. The GLM (generalized linear model) models were as follows: where y is a vector of phenotypic observations, a is a vector of SNP effects, is a covariate of population structure (PC) effects, and is residuals following Nð0; r 2 IÞ. S and Q are incidence matrices of 1s and 0s relating y to a; , respectively. The MLMs (mixed effects linear models) were fitted as follows: where y is a vector of observed phenotype, b is a vector containing coefficients of the fixed effects (SNP effects and population structure effects), u is a vector of random additive effect for individuals with Var u ð Þ ¼ r 2 g K, where K is the kinship matrix. X and Z are incidence matrices mapping y to b and u. e is the random residual effect following Nð0; r 2 e IÞ (Bradbury et al. 2007).

Identification of Top Candidate Genes
Because natural selection is likely to lead to elevated LD across relevant genes or regulatory regions, we performed a gene-based analysis of SNP-phenotype associations to identify candidate genes or genomic regions underlying adaptive traits for both transects. We first annotated all SNPs to genes that contained them or to genes within 2 kb based on the P. trichocarpa v3.0 genome annotation. We then collapsed intergenic regions into 5-kb clusters to capture SNPs that fall outside the neighborhood of any gene. This resulted in binning of the 1.3 million SNPs into 42,970 genes or intergenic regions. For each trait, we identified SNPs in the first percentile of P values from the selected model as outliers. The number of SNP outliers and total number of SNPs were counted for each gene or intergenic region, and classified as top candidate genes when their proportions of significant SNPs exceeded the 0.999 quantile of the binomial expectation (Yeaman et al. 2016). The expected proportion of outliers per gene was defined as the mean proportion of SNP outliers across genes with at least five SNPs in total and with at least one SNP outlier. Finally, we assigned gene-wise P values for each gene or intergenic region based on a binomial distribution (Yeaman et al. 2016). This gene-based method has more power to detect genomic clusters of loci possessing pronounced association signals that exceed the genomic background and reduces false positives caused by isolated outliers (Liu et al. 2010;Kang et al. 2013).

Test of Parallel Evolution across Two Transects
The null-W test of convergence, first proposed by Yeaman et al. (2016), is a more powerful test for uncovering signatures of convergence, which may not be detected by comparing direct overlap in gene outliers. The null-W test adjusts for factors such as gene size and linkage that may influence the test's significance by comparing candidate and noncandidate genes to a SNP control panel (Yeaman et al. 2016). For each candidate gene in one transect, this test evaluates the association strength of the same gene in the other transect by comparing it to a null distribution. To construct the null distribution, we extracted the R 2 of a random sample of 10,000 SNPs from noncandidate genes (i.e., those not identified as "top candidates") as the control panel, and the R 2 of all SNPs annotated to noncandidate genes. We then calculated the test statistic W with a Wilcoxon-signed rank test by comparing SNP R 2 's of each noncandidate genes versus that of the SNP control panel, using the R function "wilcox.test." These 10,000 null-W test statistics were converted into a null distribution of Z scores using Z ¼ 2W À n 1 n 2 ð Þ = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n 1 n 2 ðn 1 þ n 2 þ 1Þ=3 p , where n 1 and n 2 are the number of SNPs in the SNP control panel and noncandidate genes (Whitlock and Schluter 2009;Yeaman et al. 2016). Finally, the R 2 of SNPs annotated to candidate genes were compared with the 10k control panel their W test statistic calculated and converted into Z scores. The P value for the null-W test was calculated as the number of genes with Z score exceeding the observed Z score in the null distribution divided by the total number of genes in the null distribution. The empirical P values were further adjusted using the false discovery rate (Benjamini and Hochberg 1995) and candidate genes with adjusted P value 0.05 were considered as parallel outliers.

Genotype-Environment Association Analysis
Climate variables were estimated for each sampling site with ClimateWNA software (Wang, Hamann, et al. 2012 . We used Bayenv2 (Gunther and Coop 2013) to test for univariate associations between allele frequency and the 21 environment variables, as well as latitude, altitude, and elevation. We first estimated the population covariance matrix based on a set of 10,000 intergenic SNPs, and then tested for correlations between allele frequencies and the climatic variables using 100,000 iterations, averaged across three runs. For each climate variable, we classified SNPs with average log 10 (BF) > 2 as selection outliers. We also performed a gene-based analysis as above to identify candidate genes with excessive number of climate-associated SNPs, and null-W tests for parallel adaptation were performed on top candidate genes identified in each transect. The squared spearman correlation coefficients were used in constructing the null distribution and the null-W test.
We also used redundancy analysis (RDA), a multivariate constrained ordination method, to test for gene-environment relationships (Forester et al. 2016(Forester et al. , 2018. RDA fits multivariate linear regressions between genotype and environment to obtain a matrix of fitted values, and then uses principal component analysis on the fitted values to produce ordination axes. To reduce the multicollinearity among climate variables, we estimated correlations among them and retained seven representative and relatively uncorrelated variables (MAT, MWMT, TD, MAP, AHM, SHM, FFP, PAS, EMT, EXT, and Eref). For each of the first three axes, SNPs with a "locus score" exceeding the mean 6 3 SD were considered outliers associated with their most highly correlated climate variable. All RDA calculations were performed with the R package "vegan" (Forester et al. 2018).

Annotation of Top Candidate Genes
To identify overrepresented biological processes or molecular functions among top candidate genes identified for each trait, we first performed gene ontology (GO) enrichment analysis with R package topGO. Genes with at least five SNPs and at least one SNP outliers were used as background sets for each trait and each transect. The parent-child method is employed in the GO enrichment analysis to account for the dependency structure among GO terms. Only GO terms associated with at least five genes were tested for GO enrichment and those with P value < 0.05 were considered as interesting biological processes. The same GO enrichment analysis procedure was performed on parallel adaptation outliers with all candidate genes as the background gene set.

Clinal Variation of Phenotypic Traits and Climatic Variables
We

Parallel Phenotype-Genotype Associations between Spatial Transects
Using the null-W test, between 3% and 10% of candidate genes showed signatures of parallel association for both transects (table 1). Genes associated with timing of bud set and bud flush had the most significant overlap across two transect (P < 1.0e À4 in all hypergeometric tests). Among shared gene outliers, Potri.016G088300 (NAM; No apical meristem) and Potri.016G088800 (transducin family protein) were associated with timing of bud break in both transects, whereas Potri.001G464700 (FAD-binding Berberine family protein) was associated with plant height and diameter in both transects (supplementary table S6, Supplementary Material online).

Genotype-Environment Associations across Two Transects
Bayenv2 analysis suggested a number of candidate genes or genomic regions strongly correlated with one or more of the 21 climate variables and 3 geographic variables ( fig. 6; supplementary figs. S10 and S11, Supplementary Material online; table 1) for either the latitude or the altitude transect. Among genes associated with mean annual temperature (MAT), GO terms including cellular response to stress and carbohydrate metabolic process were enriched across altitude, and steroid biosynthesis and response to cold across latitude (supplementary tables S7 and S8, Supplementary Material online). More than 40% of candidates associated with MWMT and MCMT were also associated with MAT across both transects. Candidate genes associated with these temperature variables across altitude included Potri.001G261300 (LRR and NB-ARC domains-containing disease resistance protein), Potri. 001G444700 (NB-ARC domain-containing disease resistance protein), and Potri.015G110100 (FRI; FRIGIDA-like protein).
Candidate genes associated with mean annual precipitation (MAP) were mostly involved in cellular component organization, histone modification, and establishment of protein localization (table 1; supplementary tables S7 and S8, Supplementary Material online). The MAP-associated genes had <4% overlap with temperature-related climate variables across both transects. Heat-moisture related genes were enriched in response to abiotic stimulus and polysaccharide biosynthetic process across altitude, and cellular aldehyde metabolic process and response to endogenous stimulus across latitude (supplementary tables S7 and S8, Supplementary Material online). For altitude, one SHMassociated candidate gene-Potri.014G170600 (COL9; CONSTANS-like 9)-may downregulate expression of CONSTANS and FT and control bud dormancy in poplar (Bö hlenius et al. 2006). In general, pathways including cellular response to stress (GO:0033554), lipid binding (GO:0008289), response to cold (GO:0009409) were overrepresented among candidate genes associated with most climate variables (supplementary tables S7 and S8, Supplementary Material online).

Parallelism in Climate-Selected Loci between Transects
We detected a significant number of candidates that directly overlapped based on the binomial test in each transect. These were associated with MCMT, TD, MAP, DD_0, and eFFP (table 1; P < 1.0e À4 based on a hypergeometric test). Additionally, the Null-W test revealed additional candidate . Bud phenology refers to candidate genes associated with timing of bud set and/or bud flush. Regeneration refers to candidate genes associated with regeneration height and/or regenerated branch number. Height-diameter refers to candidate genes associated with height and/or diameter. genes with signatures of climatic selection in both transects (supplementary tables S9 and S10, Supplementary Material online). For example, nine outliers in the altitude transect were significant based on the null-W test in the latitude transect, including Potri.011G136100 (DAWDLE; DDL), which was associated with AHM, MSP, PAS, and SHM. This gene encodes a protein similar to DAWDLE, and its ortholog in Arabidopsis is involved in biogenesis of miRNAs (Zhang et al. 2018). Among climate-associated genes identified in the latitude transect, 27 displayed parallel association with the same climate variable across altitude. One gene Potri.001G282300 (UGT87A2; UDP-GLUCOSYL TRANSFERASE 87A2) associated with mean monthly temperature difference (TD), encodes a putative glycosyltransferase, which regulates abiotic stress response and flowering time via FLC (FLOWERING LOCUS C) in Arabidopsis (Wang, Jin, et al. 2012;Li et al. 2017).

Multilocus Genotype-Climate Association
Using RDA, we identified the major climate agents driving the formation of ecotypes among latitudinal and altitudinal samples. Across the latitude transect, California and Oregon sampling locations were positively related to high MAT, high heat-moisture, and growth period (positive MAT, MWMT, AHM, SHM, and FFP); Alasak sites were characterized by low annual temperature and frost free period; northern and interior British Columbia were correlated with high precipitation as snow and continentality (positive TD, PAS, and negative SHM, AHM); and Washington and southern British Columbia were positively correlated with high MAP (fig. 7A). For the altitude transect, high elevation sites were positively related to precipitation as snow and continentality, similar to Alaska ( fig. 8A). In contrast, middle and low elevations had higher MAT, MAP, and heat-moisture.
On the first three RDA axes, we identified 16,621 candidate SNPs across altitude and 22,693 across latitude (figs. 7B and 8B). In general, SNP outliers on altitude axis one reflected associations with heat-moisture and precipitation; outliers on axis two reflected continentality, extreme temperature and frost free period; and outliers on axis three reflected temperature (MWMT, MAT) and precipitation as snow ( fig. 8B and  supplementary fig. S12, Supplementary Material online). Across latitude, candidate SNPs on axis one reflected associations with heat-moisture, MWMT, and Eref; SNPs on axis two were associated with temperature and precipitation; and SNPs on axis three were associated with TD and FFP ( fig. 7B and supplementary fig. S12, Supplementary Material online).
Among SNP outliers detected across altitude, MAT, MWMT, MAP, and TD had the highest number of candidate loci. For latitude, MAT, SHM, FFP, and Eref had the most correlated candidate SNPs. There were 2,710 overlapping SNP outliers found in both transects, among which 24 were identified based on the same environmental predictor in both transects. One SNP associated with SHM in both transect was located in an intron of Potri.002G180800 (LATE ELONGATED HYPOCOTYL; LHY), which encodes a MYB-related protein regulating circadian clock by interacting with CCA1 in Arabidopsis (Kim et al. 2003). In another case, four outlier SNPs strongly associated SHM across latitude were also correlated with MAT, FFP, EMT, and PAS across altitude. These four SNPs are located on chromosome 10 within the introns of gene Potri.010G179700, which encodes flowering regulator FT2 (FLOWERING LOCUS T). This gene was previously found regulating growth and dormancy cycling in poplar (Hsu et al. 2011;Wang et al. 2018).

Overlapping Signals between Genotype-Phenotype and Genotype-Climate Associations
There was a significant overlap between gene outliers from GWAS and Bayenv2 for both altitude (5.67% of GWAS and 8.85% of Bayenv were in common; P ¼ 4.91e À15 based on a hypergeometric test) and latitude (2.78% of GWAS and 9.73% of Bayenv outliers were in common; P ¼ 2.16e À9 based on hypergeometric test) (supplementary fig. S13, Supplementary Material online). We also observed substantial overlap between SNP outliers detected with the two genotype-environment association analysis approaches. Among candidate SNPs identified by RDA, 47.3% were in common with Bayenv2 results for altitude and 55.0% for latitude. Among the eight selected climate predictors, MAT and MAP had the most significant overlap in SNP outliers detected in both Bayenv2 and RDA across both transects (P < 2.53e À5 in all hypergeometric tests; supplementary table S11, Supplementary Material online).

Phenotypic Variation Reflects Adaptation to Climate Regimes
Our results support previous findings of adaptive phenotypic gradients among natural populations along the latitudinal species range of P. trichocarpa (Evans et al. 2014;McKown, Guy, et al. 2014), and we further demonstrated that clines related to altitudinal adaptation also drive phenotypic differentiation. The strong phenotypic differentiation among natural accessions of P. trichocarpa along altitudinal and latitudinal gradients suggests that divergent selection related to local climate shapes phenotypic variation in the seasonal growth, cold hardiness, and regeneration ability of local populations across both coarse and fine spatial scales. Among the traits we measured, timing of bud phenology is a key adaptive response to seasonal variation in temperature. Provenances from high latitudes usually require a longer critical day length to induce bud formation and thus have earlier growth cessation and timing of bud set (Hurme et al. 1997;Rohde et al. 2011). In spring, bud break requires autumn chilling between 0 and 5 C, followed by sufficient heat-sum (Campbell and Sugano 1975;Hunter and Lechowicz 1992;Li et al. 2010). The late bud flush of northern and high elevational genotypes in our study may in part be due to an unsatisfied chilling or heat-sum requirements (Li et al. 2010;Basler and Korner 2014). With shorter growing seasons after transfer to nonlocal environment, most northern and high elevation genotypes displayed less growth than those originating from the center of the species range or low elevation. Although development and release of cold hardiness depends in part on daylength, it is also controlled by temperature. The spatially varying freezing tolerance we observed also reveals adaptation to different temperature regimes. Finally, regeneration ability is an important fitness measure for species that reproduce partially through clonality (Mart ınkov a and Klime sov a 2016), but less is known about the molecular determinants of these traits. The clinal pattern of regeneration ability may reflect adaptive differentiation in regeneration potential (Ikeuchi et al. 2016).

Phenotype-Associated Genes Conserved across Geographic Transects
We uncovered numerous loci with annotations suggestive of their involvement in timing of bud phenology, height, diameter, regeneration height, and branch number. The gene-based approach we employed captured genomic regions harboring concentrations of significant SNP associations in excess of the genomic background. Compared with traditional GWAS, gene-based analysis aims to identify clusters of SNP outliers and thus is less prone to false positives. Annotations of some candidate loci suggested their possible roles in controlling the corresponding phenotypes. For example, Potri.015G107000 (Regulator of chromosome condensation [RCC1] family protein), which was associated with timing of bud set, may be involved in cold acclimation (Ji et al. 2015). Another gene associated with timing of bud flush, Potri.016G088300 (No apical meristem protein; NAM) encodes a key regulator in shoot apical meristem determination (Souer et al. 1996).
By assessing the direct overlap between candidate genes, we observed common phenotype-associated genes shared by the two transects. Among them, genes associated with timing of bud set and bud flush in Virgina garden showed the most overlap in the hypergeometric test. The null-W test results provided limited evidence of shared genetic architecture underlying phenotypic adaptation to our two geographic clines. This lack of parallelism may due to divergent selection forces targeting different genes or genomic regions between latitude and altitude. However, we cannot rule out the possibility that limited sample size and control for population structure may misclassify some true positives and reduce the power of the null-W test.

Parallel and Discordant Patterns of Climate-Related Selection
The latitude and altitude transects we studied represent similar environmental clines that vary in the type and amount of precipitation, mean and seasonal temperatures, and yearly timing of growing seasons. One major shared selection agent across these two transect is local differences in MAT. Northern genotypes and those from high altitude are adapted to lower MAT and colder winters. Consequently, we expect the genomes of local populations that experience similar temperature regimes to respond in parallel, to the extent that genetic constraints exist for the loci of adaptation. To detect these signatures, we performed Bayenv2 and RDA scans and discovered a significant level of shared outlier SNPs and genes between the two transects, which is unlikely to have occurred by chance. This suggests that local selection may target similar standing genetic variations and drive adaptive allele frequency shifts across spatial groups. The GO enrichment analysis also pointed out three key pathways-carbohydrate metabolic process, cellular component organization, and lipid binding-enriched across both transects. The climate variables that showed the greatest degree of overlap between transects were related to temperature and precipitation, whereas genotype-environment relationships for variables related to evaporative demand (i.e., the interaction between temperature and precipitation in summer) did not significantly overlap. This dichotomy likely reflects the strong gradient in evaporative demand across latitude, and the lack of such a gradient across the altitudinal transect.
Despite the parallel selection signals, we also observed a large number of discordant selection outliers between the latitude and altitude transects. This is likely in part due to selective constraints unique to each transect favoring different genetic variants. For example, atmospheric pressure and daily temperature difference tend to influence plant physiology and metabolism across altitudinal gradients, whereas day length and solar angle of incidence vary across latitudinal gradients (Altshuler and Dudley 2006). The RDA revealed that temperature-related variables and yearly temperature difference are key selective forces driving population differentiation across our altitude transect, whereas temperature, heat-moisture, and seasonal growth period are the prominent selection factors for latitude. With climate gradients of varied scale across the two transects, the association strength of causal genes and pathways differ correspondingly. At the same time, historical distributions and complex demographic processes specific to local populations likely also contribute to the lack of concordance between two transects (Renaut et al. 2014;Holliday et al. 2016). Future analysis should explore how genomic features and population history may limit parallelism.

Convergence and Discordance between Selection Scans
The gene-based analysis of structured GWAS and Bayenv2 results identified genes underlying phenotypic adaptation and climatic selection, respectively. Many climate selection signals overlapped with genotype-phenotype associations, which suggests that phenotypic adaptation during colonization of the two transects was driven by selection agents captured by climate variables. For example, Potri.008G166600 (CRY2; Cryptochrome 2) was associated with timing of bud set and annual heat-moisture (AHM) in the latitude transect ( fig. 9A), and a height-associated gene (Potri.015G002300; PSEUDO RESPONSE REGULATOR 5; PPR5) also contained SNPs strongly associated with MAT and TD ( fig. 9B). This latter gene was also associated with plant height in an independent GWAS study (McKown, Guy, et al. 2014). Although we identified many cases of convergence between GWAS and climate outliers, there was also extensive discordance between these two methods, which may be due to the two methods scanning the genome for different type of association signatures. That is, the traits we measured for GWAS likely did not capture the full breadth of phenotypes responsive to climatic variation across our sampling sites.
Although there was substantial concordance, many outliers were unique to Bayenv2 or RDA, which may reflect the different logic behind the two methods. Bayenv looks for single SNP allele frequency shift across spatial groups with sample size difference and population structure accounted (Gunther and Coop 2013), while RDA considers multilocus genetic variation across multivariate environment (Forester et al. 2016). However, the common selection outliers detected with different scanning methods provide reliable candidate targets for climate-related selection.

Comparison to Previous Local Adaptation Studies in Poplar and Other Tree Species
A moderate proportion of our phenotype-and climate-related genes were concordant with findings in two recent studies in P. trichocarapa ( 10A). The two genes shared by all three studies were Potri.002G165900, an ABAinsensitive like protein that suggests a response to ABA during fall bud set, and Potri.003G214200, a glucan synthase-like protein that may be essential for 1,3-beta glucan synthesis and formation of callose deposition. The latter is particularly interesting as experimental evidence suggests that short-day induced callose deposition seals cell wall plasmodesmata in the shoot apical meristem and disconnects growth regulator transport to the dormant bud, thus maintaining dormancy (van der Schoot and Rinne 2011; Paul et al. 2014). We also found overlap between our Bayenv2 outliers and those identified with Bayenv in Evans et al. (2014) (fig. 10B). The conserved environment-associated genes between these studies include Potri.009G058900 (similar to DNA topoisomerase), Potri.012G003300, and Potri.001G359100 (Glycogen branching enzyme Thus, each of these analyses addressed a different portion of the climate space occupied by P. trichocarpa, which may partly explain the outliers that were unique to each study. Finally, we compared our results with those of three additional temperate/boreal tree species that exhibit similar patterns of local adaptation to climate. Among candidate genes identified in European aspen (Populus tremula) by Wang et al. (2018), 19 out of 91 overlapped with candidates in our RDA in both altitude and latitude transects. Among these overlapping genes was Potri.010G179700.1, which encodes FT2, a PEBP family protein that was previously found regulate growth cessation in poplar (Hsu et al. 2011). The FT2 ortholog Potra001246g10694 in Populus tremula governs clinal variation in bud set timing through variation in its transcript abundance and timing of expression (Wang et al. 2018). In addition, we found six overlapping genes between the local adaptation candidates identified by Wang et al. and our GWAS results, but no overlap with our Bayenv2 outliers. We further evaluated level of overlap between our candidate genes and adaptive loci identified in the much more distantly related interior spruce (Picea glauca Â P. engelmanii) and lodgepole pine (Pinus contorta) (Yeaman et al. 2016) and found a relatively large number of overlapping candidates (supplementary table S12, Supplementary Material online). There was significant overlap between both pine and spruce gene outliers and our Bayenv2 candidate genes across altitude (P < 0.01 in both hypergeometric tests), and our RDA candidate genes also showed significant overlap across both transects with convergent outliers in both pine and spruce (P < 0.001 in hypergeometric tests). The common outliers included several genes with roles in response to abiotic stress or seasonal light signaling: Potri.003G044200 and Potri.006G224600, which are homologs of HAB1 (HYPERSENSITIVE TO ABA1); Potri.004G102100 and Potri. 017G112700, which are homologs of the chloroplast protein APE1 (ACCLIMATION OF PHOTOSYNTHESIS TO ENVIRONMENT); Potri.002G089000, which encodes PHYA (PHYTOCHROME A); and Potri.011G119500, which encodes SPA-related 3 (SUPPRESSOR OF PHYA1) (supplementary table S12, Supplementary Material online). The large number of common genes between these studies may reflect that each used the same gene-based strategy to discern candidate genes or contigs, which emphasizes the power of this approach to define constraints in the genotype-phenotypeenvironment map, even over deep evolutionary time.

Conclusions
Climate is a primary abiotic constraint directing adaptive change in morphology, phenology, and physiology of plants. We demonstrated that adaptive traits among P. trichocarpa populations vary as a function of the local climate at their geographic origin. By employing a gene-based analysis approach to GWAS, we discovered a number of candidate loci underlying timing of bud phenology, biomass, and regeneration ability and found significant genetic parallelism underlying phenotypic adaptation between latitudinal and altitudinal transects. We also used two complementary genotype-environment scans and detected regions of the genome that have been targeted by climate-related selection, particularly temperature and its interaction with heat, and some of these regions overlapped between our altitudinal and latitudinal transects. Taken together, these results suggest some constraint in the genetics of climatic adaptation across latitude and altitude in P. trichocparpa. At the same time, adaptation along these two clines involve a large number of unique loci, which may reflect differences in their respective environmental constraints, and possibly also the demographic details of postglacial recolonization.

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.