\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Evidence of climate-driven selection on both tree traits and trait
plasticity across the climatic range of a riparian foundation species}



\author[1]{Hillary Cooper}%
\author[1]{Rebecca Best}%
\author[1]{Lela Andrews}%
\author[1]{Jaclyn Corbin}%
\author[1]{Iris Garthwaite}%
\author[1]{Kevin Grady}%
\author[1]{Catherine Gehring}%
\author[2]{Kevin Hultine}%
\author[1]{Thomas Whitham}%
\author[1]{Gerard Allan}%
\affil[1]{Northern Arizona University}%
\affil[2]{Desert Botanical Garden}%


\vspace{-1em}



  
  \date{May 24, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Selection on quantitative traits by divergent climatic conditions can
lead to substantial trait variation across a species range. In the
context of rapidly changing environments, however, it is equally
important to understand selection on trait plasticity. To evaluate the
role of selection in driving divergences in traits and their associated
plasticity within a widespread species, we compared molecular and
quantitative trait variation in Populus fremontii (Fremont cottonwood)
populations throughout Arizona. Using SNP data and genotypes from 16
populations reciprocally planted in three common gardens, we first
performed QST-FST analyses to detect selection on traits and trait
plasticity. We then explored the mechanistic basis of selection using
trait-climate and plasticity-climate regressions. Three major findings
emerged: 1) There was significant genetic variation in traits expressed
in each of the common gardens and in the phenotypic plasticity of traits
across gardens. 2) Based on QST-FST comparisons, there was evidence of
selection in all traits measured; however, this result varied from no
effect in one garden to highly significant in another, indicating that
detection of past selection is environmentally dependent. We also found
strong evidence of divergent selection on plasticity across environments
for two traits. 3) Traits and/or their plasticity were often correlated
with population source climate (R2 up to 0.77 and 0.66, respectively).
This suggests that steep climate gradients across the Southwest have
played a major role in shaping the evolution of divergent phenotypic
responses in populations and genotypes now experiencing climate change.%
\end{abstract}%



\sloppy


\textbf{Title} : Evidence of climate-driven selection on both tree
traits and trait plasticity across the climatic range of a riparian
foundation species

\textbf{Authors:} Hillary F. Cooper\textsuperscript{1,2}, Rebecca J.
Best\textsuperscript{3}, Lela V. Andrews\textsuperscript{1}, Jaclyn P.M.
Corbin\textsuperscript{1,2}, Iris Garthwaite\textsuperscript{3}, Kevin
C. Grady\textsuperscript{4}, Catherine A. Gehring\textsuperscript{1,2},
Kevin R. Hultine\textsuperscript{5}, Thomas G.
Whitham\textsuperscript{1,2}, Gerard J. Allan\textsuperscript{1,2}

\textbf{Author Affiliation:}

\textsuperscript{1}Department of Biological Science, Northern Arizona
University, Flagstaff, AZ 86011, USA\textsuperscript{2}Center for
Adaptable Western Landscapes, Northern Arizona University, Flagstaff, AZ
86011, USA

\textsuperscript{3}School of Earth and Sustainability, Northern Arizona
University, Flagstaff, AZ 86011, USA

\textsuperscript{4}School of Forestry, Northern Arizona University,
Flagstaff, AZ 86011, USA

\textsuperscript{5}Department of Research, Conservation and Collections,
Desert Botanical Garden, Phoenix, AZ, 85008, USA

\textbf{Corresponding author:} Hillary F. Cooper; \emph{telephone} :
(510) 735-6130, \emph{email} : Hillary.Cooper@nau.edu

\textbf{Keywords:} divergent selection, phenotypic plasticity,
Q\textsubscript{ST}-F\textsubscript{ST}, \emph{Populus fremontii} ,
local adaptation, climate change\textbf{Abstract:}

Selection on quantitative traits by divergent climatic conditions can
lead to substantial trait variation across a species range. In the
context of rapidly changing environments, however, it is equally
important to understand selection on trait plasticity. To evaluate the
role of selection in driving divergences in traits and their associated
plasticity within a widespread species, we compared molecular and
quantitative trait variation in \emph{Populus fremontii} (Fremont
cottonwood) populations throughout Arizona. Using SNP data and genotypes
from 16 populations reciprocally planted in three common gardens, we
first performed Q\textsubscript{ST}-F\textsubscript{ST} analyses to
detect selection on traits and trait plasticity. We then explored the
mechanistic basis of selection using trait-climate and
plasticity-climate regressions. Three major findings emerged: 1) There
was significant genetic variation in traits expressed in each of the
common gardens and in the phenotypic plasticity of traits across
gardens. 2) Based on Q\textsubscript{ST}-F\textsubscript{ST}
comparisons, there was evidence of selection in all traits measured;
however, this result varied from no effect in one garden to highly
significant in another, indicating that detection of past selection is
environmentally dependent. We also found strong evidence of divergent
selection on plasticity across environments for two traits. 3) Traits
and/or their plasticity were often correlated with population source
climate (R\textsuperscript{2} up to 0.77 and 0.66, respectively). This
suggests that steep climate gradients across the Southwest have played a
major role in shaping the evolution of divergent phenotypic responses in
populations and genotypes now experiencing climate change.

\textbf{Introduction:}

Understanding the processes shaping phenotypic diversification in nature
is a central objective of ecology and evolutionary biology (Schluter
2000; Bolnick \emph{et al.} 2011). Trait variation within widespread
species can be extensive due to historic demographic processes and
spatially and temporally heterogeneous landscapes exerting different
selection pressures across a species' range (Whitlock 2008). Over time,
subpopulations can become genetically and phenotypically differentiated
due to neutral processes, such as drift, gene flow, and mutation, as
well as the adaptive process of natural selection (Wright 1931; Spitze
1993; Holsinger \& Weir 2009; Leinonen \emph{et al.} 2013). Natural
selection acts on both phenotypes and phenotypic plasticity, defined as
the range of phenotypes a single genotype can express as a function of
environmental change (Nicotra \emph{et al.} 2010). The strength and
direction of selection may vary, creating a mosaic of trait means and
differences in trait plasticity across species' distributions (Chevin \&
Lande 2011). Adaptive evolution of individual traits and associated
plasticity can therefore differentially affect a population's
persistence on the landscape under a changing climate, as both alter the
range of phenotypes a population can express (Kelly 2019).

Phenotypic divergence is particularly evident in long-lived forest
trees, which often show strong genetic differences and local adaptation
among populations with ecological and evolutionary consequences for
associated species and communities (Savolainen \emph{et al.} 2007;
Leimu\emph{et al.} 2008; Hereford 2009; Whitham \emph{et al.} 2020). One
common hypothesis for the origins of phenotypic variation in trees is
local adaptation in response to climate. For example, studies
on\emph{Populus} have evidence of adaptive differences among populations
in growth, phenology, and physiological traits (Frewen \emph{et
al.}2000; Fischer \emph{et al.} 2017; Blasini \emph{et al} . 2020), and
the evolution of regionally adapted ecotypes (Grady \emph{et al.} 2011;
Evans \emph{et al.} 2014; McKown \emph{et al.} 2014; Ikeda \emph{et
al.}2017; Cooper \emph{et al.} 2019; Bothwell \emph{et al} . 2021). To
definitively show that phenotypic variation among populations is due to
divergent selection by their home climate, we need approaches that
integrate molecular and phenotypic assessments of replicated genotypes
across multiple common garden environments.

In addition to understanding the role of natural selection in shaping
trait differences, there has been an increasing interest in
understanding if and how selection acts on phenotypic plasticity itself
(Josephs 2018; Arnold \emph{et al.} 2019). Phenotypic plasticity is
expected to evolve proportionally to the variability and predictability
of the environment, with higher plasticity correlated with more
predictable and more heterogeneous environments (Lande 2009;
Lind\emph{et al.} 2011). For example, De Kort \emph{et al.} (2020) found
plastic responses to drought in woodland strawberry were higher in
topographically variable sites, while Leung \emph{et al.} (2020)
experimentally determined that plasticity evolved to a lower degree in
populations of a microalga experiencing less predictable salinity
conditions after 500 generations. The evolution of decreased plasticity
in homogeneous environments could occur when there is a net cost to
maintaining plasticity (DeWitt \emph{et al.} 1998). Plasticity is also
thought to increase in populations adapted to more benign climates
relative to harsh ones because the fitness cost of maladaptive
plasticity producing phenotype-environment mismatches will be greater
when resource availability is limited (Albert \& Simms 2002). This has
been demonstrated in studies where lower elevation plants produced
stronger plastic responses to drought compared to plants from harsher,
high-elevation sites (Gugger \emph{et al.} 2015; Akman \emph{et
al.}2021). Higher plasticity under milder (yet variable) conditions may
occur under the normal range of background environmental fluctuations,
however when extreme events occur outside of this range, theory predicts
rapid evolution of plasticity (Lande 2009). Evolution of increased
plasticity after extreme environmental shifts can allow mean phenotypes
to approach new optima by accelerating phenotypic adaptation, and may
enhance population persistence (Lande 2009; Chevin \& Lande 2010).
However, this depends on the shape of the reaction norm and the genetic
variance and covariances available for selection to act upon after the
extreme event (Chevin \& Hoffman 2017). Finally, although plasticity is
commonly studied on traits in isolation, species often respond to
changes in environmental with phenotypic plasticity in multiple traits,
termed multivariate plasticity (Schlichting 1989; Nielsen \& Papaj
2022). Plasticity in one trait can therefore alter plasticity in another
trait, changing the optimal multivariate plastic response and fitness
outcome to the new environment (Nielsen \& Papaj 2022). Correlations
among trait plasticities may constrain the evolution of plasticity,
resulting in discrete phenotypic strategies or solutions (Schlichting
1989). Together, these processes could combine to generate a
heterogeneous landscape, where selection gradients can produce marked
differences in plasticity along environmental clines.

A common test for whether natural selection is the mechanism responsible
for generating phenotypic divergence among populations is to compare
Q\textsubscript{ST}, the variation in quantitative traits, to
F\textsubscript{ST}, the variation in neutral genes (Wright 1951; Lande
1992; Spitze 1993). Q\textsubscript{ST} is the quantitative genetic
analog to F\textsubscript{ST} and measures the proportion of additive
genetic variance in a trait attributed to among-population differences.
If Q\textsubscript{ST} \textgreater{} F\textsubscript{ST}, there is
evidence that directional selection is responsible for population-level
divergence with respect to a trait of interest. If Q\textsubscript{ST} [?]
F\textsubscript{ST}, the null model that population differences are due
to genetic drift alone cannot be rejected. Finally, if
Q\textsubscript{ST} \textless{} F\textsubscript{ST}, this suggests
uniform or stabilizing selection acting to constrain among-population
divergence (Spitze 1993). Selection is expected to be uniform when
populations share the same phenotypic optimum and divergent when
phenotypic optima vary, such as across heterogeneous environments (Le
Corre \& Kremer 2012). The surge in both experimental and theoretical
Q\textsubscript{ST}-F\textsubscript{ST}studies has revealed a major role
of natural selection shaping intraspecific variation in quantitative
traits (McKay \& Latta 2002; Leinonen \emph{et al.} 2008;
Leinonen\emph{et al.} 2013), with approximately 70\% of all studies
showing Q\textsubscript{ST} \textgreater{} F\textsubscript{ST}
(Leinonen\emph{et al.} 2008).

Q\textsubscript{ST}-F\textsubscript{ST} comparisons can also be used to
test for selection on phenotypic plasticity (Josephs 2018). Lind\emph{et
al.} (2011) used Q\textsubscript{ST}-F\textsubscript{ST} to test for
selection on plasticity in development time among island populations of
the common frog, \emph{Rana temporaria,} which vary in pool drying
regimes. DeKort \emph{et al.} (2016) performed a modified Bayesian
Q\textsubscript{ST}-F\textsubscript{ST} analysis (Ovaskainen\emph{et
al.} 2011) to show selection on phenological plasticity in\emph{Alnus
glutinosa} across a latitudinal gradient. Alternatively, selection on
plasticity can be assessed by regressing a genotype's plasticity against
overall fitness or a fitness proxy (Pigliucci \& Schlichting 1996;
Arnold \emph{et al.} 2019). In Fremont cottonwood, for example, higher
plasticity in bud flush is associated with higher survival when
populations experience warmer temperatures (Cooper\emph{et al.} 2019).
Either approach can complement the use of environment-trait regressions
(Whitlock 2008) to test whether trait divergence among populations is
systematically related to climatic gradients as selection pressures.

The role of selection by past climatic in shaping intraspecific
variation in foundation species is especially important to quantify in
the American Southwest, where the effects of climate change are
pronounced (Garfin \emph{et al.} 2013; Williams \emph{et al.} 2020). The
loss of foundations species, defined as species that create locally
stable conditions, can have dramatic effects on fundamental ecosystem
processes like energy fluxes and biodiversity (Ellison \emph{et
al.}2005). The Southwest is described as one of the most
``climate-challenged'' regions of North America, with warming
temperatures and increasing drought events already contributing to
massive forest mortality events (Breshears \emph{et al.} 2005). Fremont
cottonwood is especially sensitive to drought and high temperature,
particularly in combination, as evidenced by stand-level mortality at
the Bill Williams National Wildlife Refuge on the lower Colorado River
(Fig. 1). Recent studies by Hultine \emph{et al.} (2020a) and
Blasini\emph{et al.} (2020) suggest that these trees are at the edge of
their thermal tolerance where water is essential for evaporative
cooling. This mortality is associated with the megadrought that Williams
\emph{et al.}(2020) identify as being the second worst drought in the
past 1200 years in the American Southwest. Thus, current climatic
gradients will be exacerbated by ongoing climate change, leading to new
selection pressures on functional traits that may be locally adapted to
a narrower range of environmental conditions.

In this study, we use trait data from three experimental common gardens
to quantify divergence (Q\textsubscript{ST}) in both genotype means
within environments and genotype plasticities across environments for
five traits of \emph{Populus fremontii} . We then compare these
Q\textsubscript{ST} values to neutral genetic divergence
(F\textsubscript{ST}). Common gardens are necessary to ensure that
among-population variance components reflect genetic differences and are
not inflated by environmental effects (Leinonen \emph{et al.} 2013).
Reciprocal common gardens can reveal traits that vary across
environmental gradients as a result of phenotypic plasticity (Kawecki \&
Ebert 2004; Franks \emph{et al.} 2014). Plastic responses to
environmental stress or release from stress may mask or amplify
genetically determined trait differences that have emerged as a result
of divergent selection (Oke \emph{et al.} 2015). It is therefore
important to assess phenotypes in multiple growing conditions to see how
the environment can modify the degree to which we can detect evidence of
selection. Our use of multiple common gardens adds to the
Q\textsubscript{ST} literature by examining how the detection of trait
differences depends on environmental conditions (Akman \emph{et
al.}2021) and by allowing for Q\textsubscript{ST}-F\textsubscript{ST}
tests on trait plasticity across gardens.

Both the collection and garden locations span an elevation gradient of
almost 2000 m, consistent with the species' range and including a
difference of 12\selectlanguage{ngerman}°C mean annual temperature and \textgreater{} 500 mm in
mean annual precipitation across source locations and
\textasciitilde{}350 mm across gardens. The benefit of these
experimental gardens is enhanced by the development of genomic data
based on the identification of \textgreater{} 9000 single nucleotide
polymorphisms (SNPs) for all genotypes planted. SNPs are an ideal type
of marker for Q\textsubscript{ST}-F\textsubscript{ST} analyses because
their mutation rates and the effects of drift are considered to be more
similar to loci that control quantitative traits compared to other
molecular markers, such as hypervariable microsatellites (Edelaar \&
Bjorklund 2011). Thus, the only difference between quantitative trait
loci driving Q\textsubscript{ST} and the loci used in
F\textsubscript{ST} estimates should be that only the latter conform to
neutral molecular evolution (Leinonen \emph{et al.} 2013).

In order to address whether climate-driven natural selection drives
trait and trait plasticity divergences across the range of Fremont
cottonwood, we evaluated three hypotheses: 1) Genetic variation in tree
traits will be evident among populations and genotypes in each of the
three common gardens, although the magnitude of the genetic effects may
vary across environments and among traits. Likewise, populations will
differ in the magnitude of plasticity of these traits measured across
the garden environments 2) Q\textsubscript{ST} values will be
significantly higher than the neutral expectation of
F\textsubscript{ST}, suggesting divergent selection has outweighed drift
in shaping divergence in trait means and plasticity among populations.
3) Mean population phenotypes will show strong relationships with their
climate of origin, as is expected when climate is a primary selective
force. Similarly, trait plasticity will also be correlated with
population source climate. Such plasticity-climate relationships should
emerge when population origins differ not only in mean climate
conditions but also in climatic variability across seasons and years, as
is the case in the Southwest.

\textbf{Materials and Methods:}

Collection sites and common gardens

To establish the common gardens, 16 populations of \emph{Populus
fremontii} were collected throughout Arizona, encompassing the
environmental variation experienced by the Sonoran Desert ecotype, as
well as three populations located on the Mogollon Rim within the
Colorado Plateau region of northern Arizona (Fig. 2). These populations
group genetically with and have been alternatively identified as the
Mogollon Rim (Blasini \emph{et al.} 2020) or Utah High Plateau ecotype
(Ikeda \emph{et al.} 2017; Supplemental Table 1). This sampling design
does not include the third described ecotype of the Central California
Valley (Ikeda \emph{et al.} 2017). Cuttings were taken from individual
tree genotypes located over 20 m away from each other to ensure
independent genotype sampling. Clonal replicates from 12 trees per
population were planted in the summer and fall of 2014 in each of the
three common garden sites after rooting in the greenhouse for
approximately four months.

The three replicated experimental common gardens span broad elevation
and climatic gradients, resulting in extreme climatic transfers for some
populations. The northernmost garden represents the cold edge of the
species' climatic range. It is located adjacent to Canyonlands National
Park, Utah and is maintained by The Nature Conservancy's Dugout Ranch.
The middle Arizona garden is located adjacent to the Agua Fria River in
Agua Fria National Monument and is maintained by the Arizona Game and
Fish Department. The southernmost garden is in Yuma, Arizona near Mittry
Lake, and is maintained by the Bureau of Land Management. These gardens
span over a 1500 m elevation difference, a 12°C mean annual temperature
range (10.7°C in Yuma, 17.2°C in Agua Fria, and 22.8°C in Canyonlands),
and a precipitation difference of \textasciitilde{}350mm (Supplemental
Table 1). Each common garden was planted with 4,096 trees. These trees
were arranged into four replicated blocks to account for within-garden
environmental variance, with each block made up of 16 randomized
population-level plots. Each population plot had 64 trees, made up of
three to six replicates of the 12 genotypes collected for that
population. Plots were arranged in a randomized 8 x 8 grid, with trees
spaced 1.85m in each cardinal direction. The garden was designed using
population plots instead of fully randomized by genotype to assess
population-level effects on dependent community members such as
arthropods and mycorrhizae, as well as ecosystem-level traits like
carbon flux.

In order to examine the relationship between climate and traits, we
downloaded 30-year normals (1961-1990 means) for 21 abiotic climate
variables for each of the 16 provenance sites and the three common
gardens using the program ClimateWNA (Wang \emph{et al.} 2012). Because
variation in both temperature and precipitation in the Southwest are
very strongly correlated with elevation, these current climate variables
are excellent proxies for the climates that trees have experienced
during their local evolutionary histories (r \textgreater{} 0.985 for
correlations between current MAT and MAP (WorldClim 2, Fick \& Hijmans
2017) and those variables estimated from 6,000 or 22,000 years ago
(WorldClim 1.4, Hijmans \emph{et al.} 2005). To create a multivariate
climatic index representing the environmental variation found throughout
the 16 provenances, the ClimateWNA variables plus elevation, latitude,
and longitude, were combined in a principal component analysis (PCA)
using \emph{labdsv} (Roberts 2007) and \emph{vegan} (Oksanen \emph{et
al.} 2016) packages in the R statistical language (R Core Team 2014).

Trait analysis

We analyzed five traits for phenotypic differentiation: fall bud set,
spring bud flush, specific leaf area (SLA), height, and trunk basal
diameter. Phenology of bud set and bud flush were measured in the fall
of 2015 and the spring of 2016, respectively. Bud set was recorded as
the initiation of bud formation, where internode elongation had ceased
and the newly emerged, bundled leaves were clustered at the same level
on the stem, offset from the shoot axis (Frewen \emph{et al.} 2000). Bud
set was measured at 6-10 day intervals from September through December
of 2015 on three replicates of all 12 genotypes per population in each
garden. We scored trees based on the bud stage exhibited by 50\% or more
of the apical meristems. There was little within-plant variation in
apical bud development, so we felt this was a good approximation of
whole plant progression towards dormancy. Spring bud flush was recorded
as the date when all leaves had flushed from the tree, specifically when
a leaf had expanded from all buds on the branch. Bud flush was measured
every two weeks from February through the end of April in the Yuma and
Agua Fria gardens, and through the end of May in the Canyonlands garden,
where colder temperatures persist later into the spring. Bud set, bud
flush, and height were assessed the full first year of growth
(2015-2016) for every genotype in each of the common gardens.

Specific leaf area was measured using the average of three to six fully
expanded leaves that were free of or had minimal herbivore damage (if no
leaves without damage could be found). To standardize for light
availability, all leaves were collected from a single south-facing
branch, collected approximately at breast height in May and June of
2020. Due to mortality in the gardens over the five years since
planting, and the time required for sampling, SLA was only measured on
12 populations in Yuma and Agua Fria, and seven populations in
Canyonlands, with 3-7 genotypes measured within each population.
Although early SLA analysis from 2015 showed similar patterns, concern
over residual maternal and greenhouse effects cautioned against using
this first year of data. Leaves were scanned and leaf area was measured
using ImageJ software (Schneider \emph{et al} . 2012). After the area
scans, leaves were dried using silica beads and then weighed. Lastly, we
used trunk basal diameter recorded at the end of the fourth year of
growth (2018), which was measured as the diameter at root collar (DRC),
\textasciitilde{}10cm from the soil, on every live tree in the gardens.
Diameter at root collar was used instead of diameter at breast height
because it allows us to track tree growth consistently from planting up
to their current stature.

Genetic analysis

Genomic DNA was extracted from \textasciitilde{}0.2g silica-dried leaf
tissue from all 192 genotypes using the Thermo Scientific MagJET Genomic
DNA Kit (Thermo Scientific). Double-digest restriction-associated DNA
(ddRAD) libraries were prepared using 2-5ng of DNA per sample in 20\selectlanguage{greek}μ\selectlanguage{english}L
reactions following a modified Peterson \emph{et al.} (2012) protocol.
Briefly, restriction and ligation was carried out simultaneously in 20\selectlanguage{greek}μ\selectlanguage{english}L
reactions using restriction enzymes MspI and EcoRI and universal adapter
sequences for indexing PCR. Ligation products were amplified using 25
cycles of PCR. After indexing, products were checked on an agarose gel
and purified. Samples were then pooled and size selected for fragments
between 200 and 350bp using a Pippin Prep (Sage Science, Inc., Beverly,
MA). The size-selected pool was quantified by qPCR and sequenced on an
Illumina MiSeq Desktop Sequencer (Illumina, Inc. San Diego, CA) in 2x75
mode. Sequence reads were processed using a modified Stacks pipeline
(Catchen \emph{et al.} 2013; Andrews 2018). Potential chloroplast and
mitochondrial sequences were filtered from the dataset by comparing them
to other \emph{Populus} sequences downloaded from GenBank. Specifically,
we removed sequences that matched chloroplast sequences from \emph{P.
fremontii} and mitochondrial sequences from \emph{P. tremula} x \emph{P.
alba.} Parameter values for clustering were based on tests following
parameterization in Mastretta-Yanes \emph{et al.} (2015). The minimum
stack depth for each individual was three and the minimum number of
individuals per locus cluster was three. All loci were used in the
measure of F\textsubscript{ST} were found to be in Hardy-Weinberg
equilibrium. To calculate F\textsubscript{ST} and a 95\% confidence
interval around F\textsubscript{ST}, we bootstrapped population-level
pairwise F\textsubscript{ST} values 500 times using
the\emph{divPartCalc} function in the R package \emph{diversity}
(Keenan\emph{et al.} 2013).

Statistical Analyses

To investigate the within- and among-population and genotype variation
in phenotypic traits, each garden was modeled separately using linear
random effects models fit by maximum likelihood in the
\emph{lme4}package in R (R Core Team 2014; Bates \emph{et al.} 2015).
The tree traits were modeled as response variables, while population and
genotype nested within population were random effects. Garden plot was
included as a random variable to help account for within-garden
environmental variance. Statistical significance was calculated using
likelihood ratio tests for the random effects using the package
\emph{lmerTest}(Kuznetsova \emph{et al.} 2015).

To model population and genotype variation in phenotypic plasticity, we
first needed to obtain replicated estimates of plasticity for each
genotype and each trait. For each genotype with at least one tree in
each of all three gardens, we randomly assigned all available trees into
genotype triplets, with one tree from each garden. We then calculated
plasticity for each triplet as the absolute value of the maximum
difference in those three trait values. This produced a number of
estimates of plasticity equal to the lowest number of trees available
for that genotype in any garden. We repeated this random triplet
assignment 100 times to obtain a set of possible plasticity datasets.
For each dataset, we estimated the variance components necessary to
calculate Q\textsubscript{ST} using a linear random effects model as
described above, where trait plasticity was the response variable, and
population and genotype nested within population were the random
effects.

For each trait and trait plasticity, we compared the quantitative trait
variation (Q\textsubscript{ST}) with genetic variance at neutral loci
(F\textsubscript{ST}). To calculate Q\textsubscript{ST} we used the
following formula:

where \selectlanguage{greek}σ\selectlanguage{english}\textsuperscript{2}\textsubscript{P} is the additive genetic
variance among populations and \selectlanguage{greek}σ\selectlanguage{english}\textsuperscript{2}\textsubscript{G} is
the additive within-population variance (Spitze 1993; McKay \& Latta
2002), i.e., the variance among genotypes within populations. Each trait
or plasticity was analyzed using the models described above, and
population and genotype variances were extracted to calculate
Q\textsubscript{ST}. Parametric bootstrap and Bayesian estimation are
considered the best methods to obtain a precision estimate around
Q\textsubscript{ST} (O'Hara \& Meril\selectlanguage{ngerman}ä 2005). We performed parametric
bootstrapping to obtain a 95\% confidence interval for
Q\textsubscript{ST}, resampling the 16 populations with replacement 1000
times, and estimating Q\textsubscript{ST} for each bootstrapped data
set. Resampling over the highest level in a hierarchical experimental
design (here the population) is considered best practice (O'Hara \&
Merilä 2005). Variance in Q\textsubscript{ST} becomes quite large as the
number of populations decreases (\textless{} 20), especially if
populations are highly differentiated (O'Hara \& Merilä 2005; Goudet \&
Büchi 2006). Goudet \& Büchi (2006) recommend sampling many populations
relative to the number of families. Our design of 16 populations with 12
genotypes per population comes close to their recommended sampling
design of upwards of 20 populations with 10 families (O'Hara \& Merilä
2005; Goudet \& Büchi 2006). In using clonally replicated genotypes, our
estimate of \selectlanguage{greek}σ\selectlanguage{english}\textsuperscript{2}\textsubscript{G} includes both additive
and non-additive genetic effects, an approach that has been shown to
lower Q\textsubscript{ST} estimates and is thus a conservative test of
Q\textsubscript{ST} \textgreater{} F\textsubscript{ST} (Cubry \emph{et
al.} 2017). Conversely, lower Q\textsubscript{ST} estimates derived from
non-additive genetic effects contribute to a more liberal test of
convergent selection (Q\textsubscript{ST} \textless{}
F\textsubscript{ST}) (Whitlock 2008; Cubry \emph{et al.} 2017). To
determine whether Q\textsubscript{ST} was significantly different from
F\textsubscript{ST}, we compared the 95\% confidence intervals for both,
which provides much stronger inference than simply comparing
Q\textsubscript{ST} to the mean F\textsubscript{ST} value (Whitlock
2008; Leinonen et al. 2013). Broad-sense heritability
(H\textsuperscript{2}) point estimates and confidence intervals were
also calculated for each trait in each garden using the equation,
H\textsuperscript{2} =
\selectlanguage{greek}σ\selectlanguage{english}\textsuperscript{2}\textsubscript{G}/(\selectlanguage{greek}σ\selectlanguage{english}\textsuperscript{2}\textsubscript{G
+}\selectlanguage{greek}σ\selectlanguage{english}\textsuperscript{2}\textsubscript{E),} where
\selectlanguage{greek}σ\selectlanguage{english}\textsuperscript{2}\textsubscript{E} includes both the plot variance
and the error variance. Calculations of heritability for plasticity did
not include plot-level variance since plasticity was measured across
gardens.

In order to test whether phenotypes showed strong climatic
relationships, we regressed population trait means and trait
plasticities against the first principal component (PC1) from the
environmental PCA. We tested these regressions and calculated an
adjusted R\textsuperscript{2} using a linear model in R (R Core Team
2014). Here, we used a single estimate of plasticity for each genotype.
Using all available replicates for each genotype, we first calculated
the mean trait value for each genotype in each garden, and then
calculated plasticity as the maximum difference between gardens.
Systematic differences among populations seen in these trait-climate
correlations are another, stronger test for evidence of divergent
selection acting over genetic drift (Whitlock 2008).

\textbf{Results:}

Genetic and Phenotypic Variation

Our dataset of 192 genotypes analyzed with ddRAD yielded 9195 SNP loci.
The 16 Arizona populations show strong differentiation with an average
pairwise F\textsubscript{ST} = 0.175 and 95\% confidence interval of
0.144-0.205. Consistent with our first hypothesis, we found significant
within and among population variation for traits at each of the three
common gardens (Table 1, Fig. 3), with phenology traits exhibiting
higher differentiation at the population than the genotype level in all
but one case. For SLA, height, and diameter at root crown (DRC), the
relative contribution of population vs. genotype varied among gardens.
Traits measured in the hottest common garden (Yuma) exhibited stronger
population than genotype effects in four out of the five traits; for
SLA, the proportion of variance explained by genotype was higher than
the proportion explained by population, although it was very close (27\%
vs. 24\%, respectively; Table 1). This garden thus produced higher
values of Q\textsubscript{ST} (the proportion of the genetic variance
that is found among rather than within populations, see next section).

Phenotypic plasticity across the three gardens showed significant
population differences in all traits except SLA (Table 1, Fig. 3). Here,
the lower sample size of seven populations with no populations from the
hottest locations may have contributed to this non-significant effect.
Population explained more variance in plasticity in the growth traits
compared to the phenology traits (Table 1).

Q\textsubscript{ST} - F\textsubscript{ST}: Comparison of quantitative
trait differentiation to neutral genetic expectation

We found evidence of selection (Q\textsubscript{ST} [?]
F\textsubscript{ST}) driving phenotype differences in over half of the
traits measured across the three gardens (Table 2, Fig. 4). Overall, the
mean Q\textsubscript{ST} value across all traits and all gardens (0.42)
was above the F\textsubscript{ST} confidence interval (0.144-0.205),
consistent with directional selection shaping trait differences and
increasing local adaptation among these populations. Bud flush traits
exhibited the highest levels of population differentiation. For example,
Q\textsubscript{ST} for bud flush in the hot and mid gardens was 0.87
and 0.86, respectively, while Q\textsubscript{ST} values for bud set
were more moderate (0.26-0.54). The confidence intervals for
Q\textsubscript{ST} crossed those of F\textsubscript{ST} in two
phenology measurements (bud set in Agua Fria and bud flush in
Canyonlands), suggesting no difference from the neutral expectation of
genetic drift of these traits in these environments. Divergent selection
on specific leaf area was apparent in the mid (Q\textsubscript{ST} =
0.54) and cold gardens (0.79), but not detectable in the hot garden of
Yuma (0.31). Tree growth traits showed relatively lower
Q\textsubscript{ST} values compared to leaf traits across gardens. Tree
height showed significant divergent selection when measured at the
hottest common garden in Yuma, Arizona (Q\textsubscript{ST} = 0.44), but
was not statistically different from F\textsubscript{ST} in the two
cooler gardens of Agua Fria and Canyonlands (0.14 and 0.26,
respectively). Contrary to our hypothesis, we found evidence of
stabilizing selection for basal trunk diameter in the coldest garden,
where the Q\textsubscript{ST} value (0.01) fell below the
F\textsubscript{ST} confidence interval. However, this result should be
interpreted with caution since dominance reduces estimates of
Q\textsubscript{ST} and is therefore a poor indicator of stabilizing
selection (Cubry \emph{et al.} 2017). This trait was indistinguishable
from the neutral expectation of F\textsubscript{ST} in the warm and mid
gardens (Q\textsubscript{ST} = 0.50 and 0.03, respectively).

Mean Q\textsubscript{ST} values for all trait plasticities except DRC
were above the F\textsubscript{ST} confidence interval (Table 2),
suggesting overall divergent selection acting on plasticity. We found
the strongest evidence for divergent selection on plasticity for bud
flush (mean Q\textsubscript{ST} = 0.84) and height (mean
Q\textsubscript{ST} = 0.66) where the Q\textsubscript{ST} 95\%
confidence interval distribution never crossed the
F\textsubscript{ST}confidence interval (Fig. 5). In the analyses of bud
set (mean Q\textsubscript{ST} = 0.44) and SLA (mean Q\textsubscript{ST}
= 0.69), the lower Q\textsubscript{ST} confidence interval distribution
overlapped with F\textsubscript{ST}, indicating that for some of the 100
possible plasticity datasets, trait plasticity differences among
populations were not distinguishable from the neutral expectation of
drift. Specifically, 67\% of the permutations fell within the
F\textsubscript{ST} interval for bud set, while only 6\% overlapped for
SLA. Similar to DRC, the Q\textsubscript{ST} values for DRC plasticity
were much lower than the other plasticities (mean Q\textsubscript{ST} =
0.07), however the upper confidence interval did overlap with
F\textsubscript{ST} 76\% of the time (Fig. 5). These three cases of
plasticity Q\textsubscript{ST} confidence interval distributions
overlapping with F\textsubscript{ST} represent a weaker detection of
selection. However, these results do provide partial evidence for SLA,
where 94/100 permutations were non-overlapping.

Climate as an agent of selection

The strength of the correlations between traits and provenance climate
varied across gardens (Table 3, Fig. 3), supporting our third hypothesis
of strong associations between phenotype and climate for some traits in
some environments. A single axis (PC1) explained 95.8\% of the variation
in provenance climate and was influenced primarily by temperature and
growing season-related climate variables (Supplemental Table 2).
Populations sourced from areas with higher temperatures, lower
precipitation, lower elevation, and longer growing seasons had higher
PC1 scores. Bud set exhibited the strongest relationship with provenance
climate across the gardens (R\textsuperscript{2} = 0.67 to 0.77), while
bud flush showed significant correlations in the two warmer gardens
(R\textsuperscript{2} = 0.49 in Yuma and 0.66 in Agua Fria), but not in
the cold garden (R\textsuperscript{2} = -0.01) (Table 3, Fig. 3a). In
Canyonlands, population variation was constrained as all trees flushed
at approximately the same time, late in the spring. Specific leaf area
did not show significant trait-climate correlations in any garden,
although we see overall SLA values increasing from the hot to the cold
garden site (Fig. 3a).

Tree growth traits were more likely to show garden-dependent
relationships between population origin and performance (Table 3, Fig.
3b). Tree height showed no climate relationships when planted in the
hottest garden, however the correlation become stronger in the mid to
cold gardens. When planted at the coldest garden, trees sourced from
colder, wetter climates, including the three populations from the
Colorado Plateau, were taller than populations from hotter, drier
environments. Similarly, the DRC relationship in the hottest garden
showed trees sourced from warm, hot environments had larger trunk
diameters compared to trees from colder climates (Fig. 3b). Together,
tree height and basal trunk diameter act as indicators that overall tree
performance is consistent with local adaptation, with hot, southern
populations growing larger in the hottest Arizona garden, and northern,
cold adapted populations growing larger in the coldest Utah garden.

Phenotypic plasticity was significantly correlated with population
source climate for all traits except SLA (Table 3, Fig. 3). Populations
sourced from hot, dry climates exhibited increased plasticity in
leaf-level phenology traits relative to the colder populations, as
previously reported in Cooper \emph{et al.} (2019). Tree-level growth
traits showed the opposite pattern of increased plasticity in those
populations sourced from the colder, high elevation environments. SLA
showed a similar trend to phenology, with warmer populations exhibiting
higher plasticity compared to populations sourced from cooler climates,
but was not significant. Again, this may be due in part to the lower
sample size of seven populations that did not include populations from
the hottest, driest sites (Fig. 3a).

\textbf{Discussion:}

We found evidence of selection acting on both traits and trait
plasticity across our three common gardens. In addition, regressions
between traits and provenance climate indicate that much of the
selection detected with Q\textsubscript{ST}-F\textsubscript{ST} analysis
is driven by climatic clines to which these populations are locally
adapted. The result of mostly high Q\textsubscript{ST} values for traits
is consistent with a majority of studies finding Q\textsubscript{ST} is
generally larger than F\textsubscript{ST} for ecological traits (McKay
\& Latta 2002; Savoleinen \emph{et al.} 2007; Leinonen \emph{et
al.}2013). However, our result of divergent selection acting on trait
plasticity is quite striking in light of the relatively few examples of
selection on plasticity documented in the literature (Arnold \emph{et
al.} 2019). Combining Q\textsubscript{ST}-F\textsubscript{ST} analysis
on trait plasticity with plasticity-climate regressions can help uncover
the evolutionary forces shaping plasticity differences across
environmental gradients (Whitlock 2008; Josephs 2018; Kelly 2019).
Below, we discuss the evidence for climate-driven adaptive divergence in
traits and trait plasticity, local adaptation to climate, and the
potential consequences of both under current climate change.

Divergent selection shapes population trait and plasticity differences

We interpret cases when both the
Q\textsubscript{ST}-F\textsubscript{ST}analysis showed large divergences
from neutral expectation and phenotype-climate correlations were
significant as strong evidence for climate-driven selection. Cases with
only one of these tests showing population differences provide partial
evidence for climate-driven selection (Table 4). For instance, there
were four cases showing Q\textsubscript{ST} \textgreater{}
F\textsubscript{ST} but non-significant trait-climate correlations.
These inconsistencies between the two tests could be due to divergent
selection that is not related to the climatic gradients we tested. There
were also four cases showing Q\textsubscript{ST} [?] F\textsubscript{ST}
and significant trait-climate correlations. In these cases,
Q\textsubscript{ST}confidence intervals overlapping with
F\textsubscript{ST} could be due to the bootstrap sampling of genotypes
with less genetic variation compared to the full sampling design, thus
lowering the Q\textsubscript{ST} estimate. There were also more
significant plasticity-climate correlations (4 out of 5) than
significant Q\textsubscript{ST}-F\textsubscript{ST} differences for
plasticity (2 out of 5). Finally, there was a range of results across
the three gardens within the Q\textsubscript{ST}-F\textsubscript{ST}
analysis itself. Together, these tests provide a continuum of support
for selection on traits and trait plasticity, and highlight which traits
may be under the strongest selection and potentially the most important
to investigate under climate change.

We found the largest Q\textsubscript{ST} values for spring bud flush,
consistent with other studies showing high phenological divergence
across latitudinal clines (Hurme 1999; Howe \emph{et al.} 2003;
Hall\emph{et al.} 2007; Evans \emph{et al.} 2016). Spring bud flush is
highly differentiated among \emph{P. fremontii} populations, with a
difference of up to eight weeks observed in flush timing (Grady \emph{et
al.} 2015; Cooper \emph{et al.} 2019; Blasini \emph{et al.} 2020). We
also found large population differences in fall bud set timing of
\textasciitilde{}2-5 weeks across the common garden gradient, reflected
in moderate Q\textsubscript{ST} values in two out of the three gardens.
The strong population differences in phenology found here agree with
Fischer \emph{et al.} (2017), who showed leaf phenology accounted for
\textgreater{}80\% of the variation in tree and forest productivity
among Fremont cottonwood genotypes. We found larger population
differences in bud flush compared to bud set. This result is intriguing
given that spring bud flush is primarily governed by temperature, while
fall bud set is mostly cued by precise day length periods (Thomashow
2001; Howe \emph{et al.} 2003). While day length is driven by latitude
and is constant from year to year, temperature can vary. The fixed
environmental cue of day length should allow populations to become
highly locally adapted and differentiated in bud set compared to a
variable environmental cue such as temperature. However, the strength of
photoperiod-driven selection on bud set may be relaxed in all but the
highest elevation populations, where the trade-off for longer growing
seasons is selected for in areas that very rarely or never experience
killing frosts (Howe \emph{et al.} 2003). Both phenology traits also
showed strong relationships with provenance climate across the gardens,
except for bud flush in the coldest garden, where low temperatures
prevented an earlier flush in the southern, warm-adapted populations
(Fig. 3a).

Our detection of selection was dependent, in part, on the environmental
conditions of each garden. Bud flush, bud set, and SLA all exhibited
divergent selection (Q\textsubscript{ST} \textgreater{}
F\textsubscript{ST}) in two out of three gardens (Fig. 4). Tree growth
traits exhibited even larger variation in Q\textsubscript{ST} among
gardens. For example, we observed high population differentiation in
height expressed in the hottest garden (Q\textsubscript{ST} = 0.44).
When populations were planted in the moderate and cool gardens, these
population differences diminished, but became more strongly predictable
from home climate (Fig 3). Q\textsubscript{ST} estimates for trunk
diameter also decreased with decreasing garden temperature. This
variability in Q\textsubscript{ST} across gardens suggests that
phenotypes shaped by selection pressures across a species' range can be
expressed differently in different growing environments, with some
environments enhancing and others dampening population phenotypic
differences (Oke \emph{et al.} 2015; Akman \emph{et al.} 2021).
Particularly for growth traits, this may represent an interaction
between the selection pressures that have shaped existing variation
across the species range and novel selection pressures imposed in a
common garden experiment or under future climate change.

The larger population-level trait differences exhibited in the hottest
common garden for most traits (except SLA) could be driven by the
maladaptation of the cold-adapted, northern populations to the extreme
thermal conditions experienced in this hot garden. This climate transfer
from northern to southern Arizona represents an extreme warming
treatment, a scenario that may be imposed on populations under severe
heat waves with climate change (Cook \emph{et al} . 2015). Similarly,
Evans \emph{et al.} (2016) found that the relationship between
Q\textsubscript{ST} and F\textsubscript{ST} changed through time, with
tree height displaying high population differentiation
(Q\textsubscript{ST} \textgreater{} F\textsubscript{ST}) under the
growing conditions in one year but not the next. Long-term common garden
experiments can demonstrate how population differences are expressed
both across different environments and through time. Given the
intensification of extreme events and climate variability going forward
(Jentsch \emph{et al.} 2007; Ganguly \emph{et al.} 2009; Garfin \emph{et
al.} 2013; Williams \emph{et al.} 2020), these types of field trials
should be expanded to evaluate the correspondence between the degree of
existing climate adaptation and the potential for future climate
survival, either through phenotypic plasticity, selection on remaining
genetic variation, or a combination of the two (Nicotra \emph{et
al.}2010; Josephs 2018).

Our Q\textsubscript{ST}-F\textsubscript{ST} comparison revealed support
for divergent selection acting on phenotypic plasticity in bud flush and
tree height, and showed partial evidence for selection on plasticity in
the other three traits (Table 2; Fig. 5). This is in contrast to
previous studies that found no evidence of selection on trait plasticity
using Q\textsubscript{ST}-F\textsubscript{ST} type comparisons
(Lind\emph{et al.} 2011, De Kort \emph{et al.} 2016), and low overall
support for selection on plastic responses to temperature (Arnold
\emph{et al.}2019). Our results suggest that for some traits,
differences in plasticity among populations across a wide environmental
gradient are larger than expected from neutral genetics, where some
populations show minimal plasticity and others exhibit high plasticity.
Conversely we found some evidence for stabilizing selection in DRC
plasticity, indicating that the difference in the magnitude of
plasticity for this trait across our populations was smaller than
expected by F\textsubscript{ST}, however it was not below the
F\textsubscript{ST}confidence interval in all 100 plasticity
permutations. The mosaic of natural selection acting on trait plasticity
across our populations shows how plasticity itself can evolve in
response to different climates.

The mosaic of natural selection acting on trait plasticity across our
populations shows how plasticity itself can evolve in response to
different climates. We found significant plasticity-climate
relationships in phenology and growth traits, where the sign of the
correlation switched between these two types of traits (Fig. 3).
Specifically, we found trees sourced from colder environments were
significantly more plastic in height and DRC compared to the warm
provenance populations, but were not as plastic with regard to their bud
set and bud flush (Fig. 3). This is an example of a multivariate
plasticity response, where plasticity in one trait may be affecting the
plasticity in another trait (Nielsen \& Papaj 2022). The higher
plasticity in phenology traits measured in populations from hotter
provenances is counterintuitive because colder source populations
experience much more predictable fall freezing events and higher yearly
temperature variation (see TD in Supplemental Table 1), and theory
predicts plasticity will increase under predictably variable
environments (Chevin \& Lande 2010). However, our climate transfer of
southern Arizona populations to the northernmost cold garden represents
an extreme climate event (over 15\selectlanguage{ngerman}°C colder in the coldest month for the
populations from the hottest source locations, Supplemental Table 1) far
outside of some populations' normal temperature range, which can result
in large, maladaptive plastic responses (Chevin \& Hoffman 2017). The
higher phenological plasticity seen in hot-adapted populations did not
translate into increased growth or growth plasticity, likely due to
maladaptive phenological plasticity that pushed these trees outside of
the appropriate growing season window (Cooper \emph{et al.} 2019). The
high bud set plasticity of warm-adapted populations meant that these
trees did not set bud until late in the growing season, when freezing
temperatures damaged non-dormant tissues. The subsequent frost damage
translated to lower growth compared to cold-adapted trees that set buds
earlier in the season and avoided frost damage. The increased height of
the cold adapted populations in the coldest garden relative to the warm
populations produced the significant differences in height plasticity.
Therefore, our result of higher height plasticity in populations sourced
from cold locations can be partially explained by the warm populations'
maladaptive plasticity in phenology.

In comparing our results to previous findings of no divergent selection
on plasticity in other systems, it is important to consider both climate
means and variances. In this system, higher growth plasticity observed
in populations sourced from colder, high elevation locations could also
be due to adaptation to increased climate variability, compared to the
central and southern Arizona populations. Specifically, the temperature
difference between the mean warmest month and the mean coldest month was
the largest for the three populations collected on the Colorado Plateau
compared to the rest of the populations below the Mogollon Rim of the
Plateau (see TD in Supplemental Table 1). This follows the theory that
higher levels of plasticity should occur in more variable environments
(Lande 2009).

Finally, our estimates of heritability in both traits and trait
plasticity also indicate that these components of the phenotype can
evolve in response to selection, at least under some environmental
conditions. Broad-sense heritability values for the five traits were
moderate, with a mean value across all gardens of 0.21 (Table 2). Our
phenology heritability measures (H\textsuperscript{2} = 0.04-0.48 for
bud flush and H\textsuperscript{2} = 0.19-0.30 for bud set) were lower
than previously found in some \emph{Populus} studies
(H\textsuperscript{2} = 0.94 for bud flush and H\textsuperscript{2} =
0.91 for bud set in \emph{P. trichocarpa x deltoides} , Frewen \emph{et
al.} 2000). Heritability values for growth traits (H\textsuperscript{2}=
0.11-0.27 for height and 0.08-0.21 for DRC) and SLA
(H\textsuperscript{2} = 0.10-0.35) were fairly consistent with other
reported \emph{Populus} estimates (H\textsuperscript{2} = 0.03-0.42 for
height and H\textsuperscript{2} = 0.09-0.25 for diameter at breast
height in \emph{P. tremuloides} , Ding \emph{et al.} 2020;
H\textsuperscript{2} \selectlanguage{english}[?] 0.2-0.6 for SLA in \emph{P. nigra} , Guet
\emph{et al.} 2015). The range of heritability estimates for the same
trait across the three gardens highlights the environment-dependent
nature of heritability. This is especially apparent in our bud flush
results, where we found the lowest value in the cold garden
(H\textsuperscript{2}= 0.04) and the highest value in the warmest garden
(H\textsuperscript{2} = 0.48). There was also no trend toward higher or
lower heritability estimates in a particular common garden. These
results suggest that in some environments evolutionary potential is
limited but can increase as environmental conditions and associated
selection pressures change. Furthermore, these heritability increases
are not necessarily associated with a specific direction of change
(\emph{i.e.,} increasing or decreasing temperature). Broad-sense
heritability for the five trait plasticities ranged from 0.09-0.18, a
similar result to bud burst plasticity found in another riparian
deciduous tree, black alder (H\textsuperscript{2} = 0-0.129, De
Kort\emph{et al.} 2016). Our results of genetic variation in trait
plasticity combined with the evidence for selection based on
Q\textsubscript{ST}-F\textsubscript{ST} analysis and non-zero
heritability estimates show selection on the heritable components of
phenotypic plasticity may lead to evolving plasticity across the
landscape among these Arizona populations of Fremont cottonwood.

Local adaptation to climate

Whereas leaf phenology and morphology traits (bud flush, bud set, and
SLA) had the highest degree of differentiation in our study, likely due
to climate-related divergent selection, tree growth traits (height and
trunk basal diameter) were less differentiated. This result suggests
that the neutral processes of gene flow and drift may override weak
selection on growth traits (McKay \& Latta 2002), or that growth is
strongly constrained by the relative success dictated by the rest of a
plant's phenotype (Saint-Laurent \emph{et al.} 2007; Leinonen \emph{et
al} . 2013). In this case, the latter explanation is very likely.
Although phenology and growth traits showed significant regressions with
climate of origin (Fig. 3), only the growth traits reversed the sign of
that relationship across gardens. Specifically, phenology trends were
mostly constant across gardens, with warmer source populations setting
bud later and flushing earlier regardless of growing environment (except
for bud flush in Canyonlands). However, height and trunk diameter
declined as transfer distance increased (in terms of both hotter and
colder climates) for populations relative to their home sites. This
indicates local adaptation, where the highest productivity is observed
in populations whose source climate best matches that of the garden
climate. In the hottest garden (Yuma), there was a positive relationship
between trunk diameter and warmer provenance climates. In the coldest
garden (Canyonlands), the reverse was true, where trees from cooler
provenances grew significantly taller than those from the warmer sites
(Fig. 3b). Whether this higher performance of local populations is
enabled by their leaf phenology and morphology traits measured here vs.
additional plant functional traits is an important area for further
study.

\emph{Management implications} \textbf{}

When local adaptation and phenotypic differentiation in forest trees are
closely tied to variation in climate, populations may become
increasingly maladapted as climate change continues (Shaw \& Etterson
2012; Franks \emph{et al.} 2014; Aitken \& Bemmels 2015). However, the
magnitude of climate change combined with the degree of genetic
variability, heritability, and phenotypic plasticity of traits will all
interact to determine the extent of adaptation or maladaptation.
Maladaptation due to climate change is expected to be greatest in
populations from the warmest extent of their range, while populations at
the cold edge may benefit from slightly warmer temperatures (Aitken \&
Bemmels 2015). This expectation corresponds with the maladaptive
phenological plasticity we have observed in southern populations vs.
adaptive phenological plasticity in northern populations (Cooper\emph{et
al.} 2019). However, nuanced changes in temperature and precipitation
patterns will produce novel genotype-climate associations, creating more
complex climate responses compared to the poleward range shifts and
vulnerable trailing edges traditionally associated with warming
(Gougherty \emph{et al.} 2021). Although our study does not encompass
the full genetic and geographic range of Fremont cottonwood, our results
of declining performance as climate transfer distance increases suggests
that this species will likely experience maladaptation as local
conditions become more arid, especially for southern populations that
are close to their thermal tolerance (Ault\emph{et al.} 2014; see Fig.
1). Because these trees are important foundation species of riparian
systems, selecting genotypes with sufficient performance under warming
conditions is essential for the persistence of associated communities
and ecosystems.

Conclusions

By combining Q\textsubscript{ST}-F\textsubscript{ST} analyses and
climate-trait regressions, we show that both quantitative traits and
associated plasticity are targets of climate-driven selection, where
population may exhibit significant differentiation in both, one, or
neither type of analysis (Table 4). This spectrum can help determine
which traits are likely under the strongest climate-drive selection.
Specifically, phenology as well as bud flush plasticity and height
plasticity show the strongest responses to selection (Table 4). These
traits and plasticities could be most impacted by further climate
change, as they show strong signals of climate-driven selection in the
past. However, both the detection of past selection and the
differentiation in current performance differed strongly across our
experimental common gardens. Strategies for management of widespread
species like Fremont cottonwood would benefit from considering the
climatic selection pressures of source locations to anticipate their
future performance, as well as the population-specific potential for
adaptive trait plasticity under changing environmental conditions.

\textbf{Acknowledgements}

This research was supported by NSF-IGERT and NSF GK-12 Fellowships (HF
Cooper), NSF Bridging Ecology and Evolution grant DEB\selectlanguage{english}-1914433 (RJ Best,
GJ Allan, R Lindroth, TG Whitham), NSF MacroSystems grant DEB-1340852
(GJ Allan, TG Whitham, CG Gehring, \& KC Grady), NSF Macrosystems grant
DEB-134056 (KR Hultine), NSF DBI-1126840 (TG Whitham), which established
the Southwest Experimental Garden Array. We thank our agency partners
for helping to facilitate use of the common gardens: Dana Warnecke and
Kelly Wolf at Arizona Game and Fish (Agua Fria), Erica Stewart at the
Bureau of Land Management (Yuma), and Barry Bakker, Phil Adams, and
Kristen and Matt Redd at The Nature Conservancy's Canyonlands Research
Center at Dugout Ranch. We acknowledge Christopher Updike, Zachary
Ventrella, Davis Blasini, Dan Koepke, and Matthew McEttrick, along with
many volunteers for help establishing and maintaining the common
gardens. We thank Helen Bothwell for her help developing and
troubleshooting the SNP library. Lastly, thanks to Jacob Cowan, Michelle
Hockenbury, Teresa Reyes, and Michelle Bem for assistance with data
collection in the field, and the Cottonwood Ecology and Community
Genetics Lab for their constructive comments and reviews.

\textbf{References}

Aitken, S.N., \& Bemmels, J.B. (2015). Time to get moving: Assisted gene
flow of forest trees. \emph{Evolutionary Applications, 9} , 271-290.

Akman, M., Carlson, J.E., \& Latimer, A.M. (2021). Climate explains
population divergence in drought-induced plasticity of functional traits
and gene expression in a South African Protea. \emph{Molecular Ecology,
30} , 255--273.

Alberto, F.J., Aitken, S.N., Alia, R., Gonzalez-Martinez, S.C.,
Hanninen, H., Kremer, A., \ldots{} Savolainen, O. (2013). Potential for
evolutionary responses to climate change - evidence from tree
populations. \emph{Global Change Biology} , \emph{19} (6), 1645--1661.

Alpert, P., \& Simms, E.L. (2002). The relative advantages of plasticity
and fixity in different environments: When is it good for a plant to
adjust? \emph{Evolutionary Ecology} , \emph{16,} 285--297.

Andrews, L.V. (2018). akutils RADseq utility: Simplified processing of
RADseq data through Stacks.\emph{Zenodo}
.~doi.org/10.5281/zenodo.1205089

Arnold, P.A., Nicotra, A.B., \& Kruuk, L.E.B. (2019). Sparse evidence
for selection on phenotypic plasticity in response to
temperature.\emph{Philosophical Transactions of the Royal Society B} ,
\emph{374} : 20180185\emph{.}

Ault, T.R.,~Cole, J.E.,~Overpeck, J.T.,~Pederson, G.T., \&~Meko,
D.M.~(2014).~Assessing the risk of persistent drought using climate
model simulations and paleoclimate data.~\emph{Journal of Climate}
,~\emph{27} , 7529--~7549.

Bates, D., Maechler, M., Bolker, B., \& Walker, S. (2015). Fitting
linear mixed-effects models using lme4. \emph{Journal of Statistical
Software, 67} , 1-48.

Blasini, D., Koepke, D., Grady, K., Allan, G., Gehring, C., Cushman,
S.A., \ldots{} Hultine, K. (2020). Adaptive trait syndromes along
multiple economic spectra define cold and warm adapted ecotypes in a
widely distributed foundation tree species. \emph{Journal of Ecology.
2020-0859.R1}

Bolnick, D.I., Amarasekare, P., Araujo, M.S., Burger, R., Levine, J.M.,
Novak, M., \ldots{} Vasseur, D.A. (2011). Why intraspecific trait
variation matters in community ecology. \emph{Trends in Ecology and
Evolution} , \emph{26,} 183-192.

Bothwell, H.M., Evans, L.M., Hersch-Green, E.I., Woolbright, S.A.,
Allan, G.J., \& Whitham, T.G. (2021). Genetic data improves niche model
discrimination and alters the direction and magnitude of climate change
forecasts. \emph{Ecological Applications, 31,} e02254

Breshears, D.D., Cobb, N.S., Rich, P.M., Price, K.P., Allen, C.D.,
Balice, R.G., \ldots{} Meyer, C.W. (2005). Regional vegetation die-off
in response to global-change-type drought. \emph{PNAS} , 102,
15144--15148.

Chevin, L.M, \& Lande, R. (2010). When do adaptive plasticity and
genetic evolution prevent extinction of a density-regulated
population?\emph{Evolution} , \emph{64,} 1143--1150.

Chevin, L.M., \& Lande, R. (2011) Adaptation to marginal habitats by
evolution of increased phenotypic plasticity. \emph{Journal of
Evolutionary Biology, 24} , 1462--1476.

Chevin, L.M., \& Hoffmann, A.A. (2017). Evolution of phenotypic
plasticity in extreme environments. \emph{Philosophical Transactions of
the Royal Society B: Biological Sciences} , \emph{372,} 20160138.

Cook, B.I., Ault, T.R., \& Smerdon, J.E. (2015). Unprecedented 21st
century drought risk in the American Southwest and Central
Plains.\emph{Science Advances} 1: e1400082.

Cooper, H.F., Grady, K.C., Cowan, J.A., Best, R.J., Allan, G.J., \&
Whitham, T.G. (2019). Genotypic variation in phenological plasticity:
Reciprocal common gardens reveal adaptive responses to warmer springs
but not to fall frost. \emph{Global Change Biology} , \emph{25,}
187-200.

Cubry, P., Scotti, I., Oddou-Muratorio, S., \& Lefevre, F. (2017).
Generalization of the Q\textsubscript{ST} framework in hierarchically
structured populations: Impacts of inbreeding and dominance.

Cushman, S.A., Max, T., Meneses, N., Evans, L.M., Ferrier, S., Honchak,
B., \ldots{} Allan, G.J. (2014). Landscape genetic connectivity in a
riparian foundation tree is jointly driven by climatic gradients and
river networks. \emph{Ecological Applications, 24} , 1000--1014.

De Kort, H. Vander Mijnsbrugge, K., Vandepitte, K., Mergeay, J.,
Ovaskainen, O., \& Honnay, O. (2016). Evolution, plasticity and evolving
plasticity of phenology in the tree species \emph{Alnus glutinosa.
Journal of Evolutionary Biology, 29,} 253-264.

De Kort, H., Panis, B., Helsen, K., Douzet, R., Janssens, S.B., \&
Honnay, O. (2020). Pre-adaptation to climate change through
topography-driven phenotypic plasticity. \emph{Journal of Ecology,
108,}1645-1474.

DeWitt, T.J., Sih, A., \& Wilson, D.S. (1998). Costs and limits of
phenotypic plasticity. \emph{Trends Ecol. Evol., 13,} 77--81.

Des Roches, S., Post, D.M., Turley, N.E., Bailey, J.K., Hendry, A.P.,
Kinnison, M.T., \ldots{} Palkovacs, E.P. (2017). The ecological
importance of intraspecific variation. \emph{Nature Ecology \&
Evolution} , \emph{2,} 57-64.

Ding C., Hamann, A., Yang, R-C., \& Brouard, J.S. (2020). Genetic
parameters of growth and adaptive traits in aspen (\emph{Populus
tremuloides} ): Implications for tree breeding in a warming
world.\emph{PLoS ONE} , \emph{15,} e022922.

Edelaar, P., \& Bjorklund, M. (2011). If F\textsubscript{ST} does not
measure neutral genetic differentiation, then comparing it with
Q\textsubscript{ST} is misleading. Or is it? \emph{Molecular Ecology}
,\emph{20,} 1805--1812.

Ellison, A.M., Bank, M.S., Clinton, B.D., Colburn, E.A., Elliott, K.,
Ford, C.R., \ldots{} Webster, J.R. (2005). Loss of foundation species:
consequences for the structure and dynamics of forested
ecosystems.\emph{Frontiers in Ecology and the Environment, 3,} 479--486.

Evans, L.M., Allan, G.J., Meneses, N., Max, T.L., \& Whitham, T.G.
(2013). Herbivore host-associated genetic differentiation depends on the
scale of plant genetic variation examined. \emph{EvolutionaryEcology,
27,} 65--81.

Evans, L.M., Slavov, G.T., Rodgers-Melnick, E., Martin, J., Ranjan, P.,
Muchero, W., \ldots{} DiFazio, S.P. (2014). Population genomics
of\emph{Populus trichocarpa} identifies signatures of selection and
adaptive trait associations. \emph{Nature Genetics} ,
\emph{46,}1089--1096.

Evans, L.M., Allan, G.J., DiFazio, S.P., Slavov, G.T., Wilder, J.A.,
Floate, K.D., \ldots{}Whitham, T.G. (2015). Geographical barriers and
climate influence demographic history and natural selection in
narrowleaf cottonwoods. \emph{Heredity, 114,} 387--396.

Evans, L.M., Kaluthota, S., Pearce, D.W., Allan, G.J., Floate, K., Rood,
S.B., \& Whitham, T.G. (2016). Bud phenology and growth are subject to
divergent selection across a latitudinal gradient in \emph{Populus
angustifolia} and impact adaptation across the distributional range and
associated arthropods. \emph{Ecology and Evolution, 6,} 4565--4581.

Fick, S.E., \& Hijmans, R.J. (2017). WorldClim 2: new 1-km spatial
resolution climate surfaces for global land areas. \emph{International
Journal of Climatology, 37,} 4302-4315.

Fischer, D.G., Wimp, G.M., Hersch-Green, E., Bangert, R.K., Leroy, C.J.,
Bailey, J.K., \ldots{} Whitham, T.G. (2017). Tree genetics strongly
affect forest productivity, but intraspecific diversity-productivity
relationships do not. \emph{Functional Ecology} , \emph{31} , 520--529.

Franks, S.J., Weber, J.J., \& Aitken, S.N. (2014). Evolutionary and
plastic responses to climate change in terrestrial plant
populations.\emph{Evolutionary Applications} , \emph{\textbf{7},}
123--139.

Frewen, B.E., Chen, T.H.H., Howe, G.T., Davis, J., Rohde, A., Boerjan,
W., \& Bradshaw, H.D. (2000). Quantitative trait loci and candidate gene
mapping of bud set and bud flush in \emph{Populus} . \emph{Genetics,154}
\textbf{,} 837--845.

Ganguly, A.R., Steinhaeuser, K., Erickson, D.J., Branstetter, M.,
Parish, E.S., Singh, N., \ldots{} Buja, L. (2009). Higher trends but
larger uncertainty and geographic variability in 21st century
temperature and heat waves. \emph{Proceedings of the National Academy of
Sciences of the United States of America} , \emph{106} , 15555--15559.

Garfin, G., Jardine, A., Merideth, R., Black, M., \& LeRoy, S.
(2013).\emph{Assessment of Climate Change in the Southwest United
States: A Report Prepared for the National Climatic Assessment} . A
report by the Southwest Climate Alliance. Washington, DC: Island Press.

Goudet, J., \& Buchi, L. (2006). The effects of dominance, regular
inbreeding and sampling design on Q\textsubscript{ST}, and estimator of
population differentiation for quantitative traits. \emph{Genetics}
,\emph{172,} 1337--1347.

Gougherty, A.V., Keller, S. R., \& Fitzpatrick, M.C. (2021).
Maladaptation, migration and extirpation fuel climate change risk in a
forest tree species. \emph{Nature Climate Change, 11} , 166-171.

Grady, K.C., Ferrier, S.M., Kolb, T.E., Hart, S.C., Allan, G.J., \&
Whitham, T.G. (2011). Genetic variation in productivity of foundation
riparian species at the edge of their distribution: Implications for
restoration and assisted migration in a warming climate. \emph{Global
Change Biology} , \emph{17} , 3724--3735.

Grady, K.C., Laughlin, D.C., Ferrier, S.M., Kolb, T.E., Hart, S.C.,
Allan, G.J., \& Whitham, T.G. (2013). Conservative leaf economic traits
correlate with fast growth of genotypes of a foundation riparian species
near the thermal maximum extent of its geographic range.\emph{Functional
Ecology} , \emph{27,} 428--438.

Grady, K.C., Kolb, T.E., Ikeda, D.H., \& Whitham, T.G. (2015). A bridge
too far: Cold and pathogen constraints to assisted migration of riparian
forests. \emph{Restoration Ecology} , \emph{23} , 811--820.

Guet, J., Fabbrini, F., Fichot, R., Sabatti, M., Bastien, C., \&
Brignolas, F. (2015). Genetic variation for leaf morphology, leaf
structure and leaf carbon isotope discrimination in European populations
of black poplar (\emph{Populus nigra} L.). \emph{Tree Physiology,
35,}850-863.

Gugger, S., Kesselring, H., Stocklin, J., \& Hamann, E., (2015). Lower
plasticity exhibited by high- versus mid-elevation species in their
phenological responses to manipulated temperature and
drought.\emph{Annals of Botany, 116,} 953-962.

Hall, D., Luquez, V., Garcia, V.M., St Onge, K.R., Jansson, S., \&
Ingvarsson, P.K. (2007). Adaptive population differentiation in
phenology across a latitudinal gradient in European aspen (\emph{Populus
tremula} , L.): A comparison of neutral markers, candidate genes and
phenotypic traits. \emph{Evolution} , \emph{61} , 2849--2860.

Hangartner, S., Laurila, A., \& Rasanen, K. (2012). Adaptive divergence
in moor frog (\emph{Rana arvalis} ) populations along an acidification
gradient: Inferences from
Q\textsubscript{ST}-F\textsubscript{ST}correlations. \emph{Evolution,
66,} 867--881.

Hereford, J. (2009). A quantitative survey of local adaptation and
fitness trade-offs. \emph{The American Naturalist} , \emph{173} ,
579--588.

Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., \& Jarvis, A.
(2005). Very high resolution interpolated climate surfaces for global
land areas.~\emph{International Journal of Climatology, 25,} 1965-1978.

Holsinger, K.E., \& Weir, B.S. (2009). Genetics in geographically
structured populations: Defining, estimating and interpreting
F\textsubscript{ST}. \emph{Nature Reviews Genetics, 10} , 639--650.

Howe, G.T., Aitken, S.N., Neale, D.B., Jermstad, K.D., Wheeler, N.C., \&
Chen, T.H.H. (2003). From genotype to phenotype: Unraveling the
complexities of cold adaptation in forest trees. \emph{Canadian Journal
of Botany} , \emph{1266,} 1247--1266.

Hultine, K.R., Allan, G.J., Blasini, D., Bothwell, H.M., Cadmus, A.,
Cooper, H.F., \ldots{} Whitham, T.G. (2020a). Adaptive capacity in the
foundation tree species \emph{Populus fremontii} : Implications for
resilience to climate change and non-native species invasion in the
American Southwest. \emph{Conservation Physiology, 8} (1),
doi.org/10.1093/conphys/coaa061

Hultine, K.R., Froend, R., Blasini, D., Bush, S.E., Karlinski, M.,
Koepke, D.F. (2020b). Hydraulic traits that buffer deep-rooted plants
from changes in hydrology and climate. \emph{Hydrological Processes}
,\emph{34} , 209-222.

Hurme, P. (1999). \emph{Genetic Basis of Adaptation: Bud Set Date and
Frost Hardiness Variation in Scots Pine} , University of Oulu Press 19.

Ikeda, D.H., Max, T.L., Allan, G.J., Lau, M.K., Shuster, S.M., \&
Whitham, T.G. (2017). Genetically informed ecological niche models
improve climate change predictions. \emph{Global Change Biology,23} ,
164--176.

Jentsch, A., Kreyling, J., \& Beierkuhnlein, C. (2007). A new generation
of climate change experiments: Events, not trends. \emph{Frontiers in
Ecology and the Environment} , \emph{5} , 365--374.

Josephs, E.B. (2018). Determining the evolutionary forces shaping G x
E.\emph{New Phytologist, 219,} 31-36.

Keenan, K., McGinnity, P., Cross, T.F., Crozier, W.W., \& Prodohl, P.A.
(2013). diveRsity: An R package for the estimation of population
genetics parameters and their associated errors. \emph{Methods in
Ecology and Evolution} , doi:10.1111/2041-210X.12067

Kelly, M. (2019). Adaptation to climate change through genetic
accommodation and assimilation of plastic phenotypes.\emph{Philosophical
Transactions of the Royal Society B: Biological Sciences} , 374,
20180176.

Kuznetsova, A., Brockhoff, P.B., \& Christensen, R.H.B.
(2015).~lmerTest: tests in linear mixed effects models. R package
version 2.0-20. https://cran. rproject. org/web/packages/lmerTest.

Lande, R. (2009). Adaptation to an extraordinary environment by
evolution of phenotypic plasticity and genetic
assimilation.\emph{Journal of Evolutionary Biology} , \emph{22} ,
1435--1446.

Lande, R. (1992). Neutral theory of quantitative genetic variance in an
island model with local extinction and colonization. \emph{Evolution,46}
, 381-389.

Le Corre, V., \& Kremer, A. (2012). The genetic differentiation at
quantitative loci under local adaptation. \emph{Molecular Ecology,
21,}1548-1566.

Leimu, R., \& Fischer, M. (2008). A meta-analysis of local adaptation in
plants. \emph{PloS one} , \emph{3} , e4010.

Leinonen, T., O'Hara, R.B., Cano, J.M., \& Merila, J. (2008).
Comparative studies of quantitative trait and neutral marker divergence:
A meta-analysis. \emph{Journal of Evolutionary Biology} ,
\emph{21,}1--17.

Leinonen, T., McCairns, R.J.S., O'Hara, R.B., \& Merila, J. (2013).
Q\textsubscript{ST}-F\textsubscript{ST} comparisons: Evolutionary and
ecological insights from genomic heterogeneity. \emph{Nature Reviews
Genetics} , \emph{14} , 179--90.

Lind, M.I., Ingvarsson, P.K., Johansson, H., Hall, D., \& Johansson, F.
(2011). Gene flow and selection of phenotypic plasticity in an island
system of \emph{Rana temporaria} . \emph{Evolution, 65} , 684-697.

Lynch, M., \& Hill, W.G. (1986). Phenotypic evolution by neutral
mutation. \emph{Evolution,} \emph{40,} 915-935.

Mastretta-Yanes, A., Arrigo, N., Alvarez, N., Jorgensen, T.H., Pinero,
D., Emerson, B.C. (2015). Restriction site-associated DNA sequencing,
genotyping error estimation and de novo assembly optimization for
population genetic inference. \emph{Molecular Ecology Resources}
,\emph{15,} 28--41.

McKay, J.K., \& Latta, R.G. (2002). Adaptive population divergence:
Markers, QTL and traits. \emph{Trends in Ecology and Evolution}
,\emph{17,} 285--291.

McKown, A.D., Guy, R.D., Klap\selectlanguage{polish}š\selectlanguage{english}te, J., Geraldes, A., Friedmann, M.,
Cronk, Q.B.C., \ldots{} Douglas, C.J. (2014). Geographical and
environmental gradients shape phenotypic trait variation and genetic
structure in \emph{Populus trichocarpa} . \emph{New Phytologist,201,}
1263--1276.

Nielsen, M.E., \& Papaj, D.R. (2022). Why study plasticity in multiple
traits? New hypotheses for how phenotypically plastic traits interact
during development and selection. \emph{Evolution} , 1--12.
doi:10.1111/evo.14464

Nicotra, A.B., Atkin, O.K., Bonser, S.P., Davidson, A.M., Finnegan, E.,
Mathesius, U., \ldots{}van Kluenen, M. (2010). Plant phenotypic
plasticity in a changing climate. \emph{Trends in Plant Science,15,}
684--692.

O'Hara, R.B., \& Meril\selectlanguage{ngerman}ä, J. (2005). Bias and precision in
Q\textsubscript{ST} estimates: Problems and some
solutions.\emph{Genetics} , \emph{171} , 1331--1339.

Oke, K.B., Bukhari, M., Kaeuffer, R., Rolshausen, G., Bolnick, D.I.,
\ldots{} Hendry, A.P. (2015). Does plasticity enhance or dampen
phenotypic parallelism? A test with three lake-stream stickleback
pairs.\emph{Journal of Evolutionary Biology, 29,} 126-143.

Oksanen, J., Blanchet, F.G., Kindt, R., Legendre, P., Minchin, P.R.,
O'Hara, R.B., \ldots{} Wagner, H. (2016). vegan: Community Ecology
Package. R package version 2.3-3.
http://CRAN.R-project.org/package=vegan

O'Neill, G.A., Hamann, A., \& Wang, T. (2008). Accounting for population
variation improves estimates of the impact of climate change on species'
growth and distribution. \emph{Journal of Applied Ecology} , \emph{45} ,
1040--1049.

O'Neill, G.A., Wang, T., Ukrainetz, N., Charleson, L., McAuley, L.,
Yanchuk, A., \& Zedel, S. (2017). \emph{A Proposed Climate-Based Seed
Transfer System for British Columbia} . Prov. B.C., Victoria, B.C.

Ovaskainen, O., Karhunen, M., Zheng, C., Cano Arias, J.M., \& Merilä, J.
(2011). A new method to uncover signatures of divergent and stabilizing
selection in quantitative traits. \emph{Genetics, 189,} 621--632

Parmesan, C. (2006). Ecological and evolutionary responses to recent
climate change. \emph{Annual Review of Ecology Evolution and Syst}
e\emph{matics} , \emph{37} , 637-- 669.

Perez-Harguindegay, N., Díaz, S.A., Garnier, E.B., Lavorel, S.C.,
Poorter, H.D., Jaureguiberry P, \ldots{} Cornelissen, J.H.C. (2013). New
handbook for standardised measurement of plant functional traits
worldwide. \emph{Australian Journal of Botany} , \emph{61,} 167-234.

Pigliucci, M., \& Schlichting, C.D. (1996). Reaction norms of
Arabidopsis. IV. Relationships between plasticity and
fitness.\emph{Heredity, 76,} 427--436.

R Core Team (2014). R: A language and environment for statistical
computing. http://www.R--project.org/

Saint-Laurent, R., Legault, M., \& Bernatchez, L. (2003). Divergent
selection maintains adaptive differentiation despite high gene flow
between sympatric rainbow smelt ecotypes (\emph{Osmerus
mordax}Mitchill). \emph{Molecular Ecology, 12,} 315--330.

Savolainen, O., Pyhäjärvi, T., \& Knürr, T. (2007). Gene Flow and Local
Adaptation in Trees. \emph{Annual Review of Ecology, Evolution, and
Systematics} , \emph{38,} 595--619.

Shaw, R.G., \& Etterson, J.R. (2012). Rapid climate change and the rate
of adaptation: insight from experimental quantitative genetics.\emph{New
Phytologist} , \emph{195,} 752--765.

Schlichting, C.D. (1989). Phenotypic Integration and Environmental
Change. \emph{BioScience} , \emph{39} , 460--464.

Schluter, D. (2000). \emph{The Ecology of Adaptive Radiation} . Oxford
University Press, Oxford.

Schneider, C.A., Rasband, W.S., \& Eliceiri, K.W. (2012). NIH Image to
ImageJ: 25 years of image analysis,~\emph{Nature Methods,~9} \textbf{,}
671-675.

Smith, D.S., Bailey, J.K., Shuster, S.M., \& Whitham, T.G. (2011). A
geographic mosaic of trophic interactions and selection: trees, aphids
and birds. \emph{Journal of Evolutionary Biology,} \emph{24} , 422--429.

Spitze, K. (1993). Population Structure in \emph{Daphnia obtusa} :
Quantitative genetic and allozymic variation. \emph{Genetics,
135,}367--374.

Thomashow, M.F. (2001). So what's new in the field of plant cold
acclimation? Lots! \emph{Plant Physiol., 125,} 89--93.

Wang, T., O'Neill, G.A., \& Aitken, S.N. (2010). Integrating
environmental and genetic effects to predict responses of tree
populations to climate. \emph{Ecological Applications} ,
\emph{20,}153--163.

Wang, T., Hamann, A., Spittlehouse, D.L., \& Murdock, T.Q. (2012).
ClimateWNA --- High-Resolution Spatial Climate Data for Western North
America. \emph{Journal of Applied Meteorology and Climatology}
,\emph{51,} 16--29.

Whitham, T.G., Bailey, J.K., Schweitzer, J.A., Shuster, S.M., Bangert,
R.K., LeRoy, C.J., \ldots{} Wooley, S.C. (2006). A framework for
community and ecosystem genetics: From genes to ecosystems. \emph{Nature
Review Genetics} , \emph{7,} 510--523.

Whitham, T.G., Allan, G.J., Cooper, H.F., \& Shuster,~S.M. (2020).
Intraspecific genetic variation and species interactions contribute to
community evolution. \emph{Annual Review of Ecology, Evolution, and
Systematics} , \emph{51} , 587-612.

Whitlock, M.C. (2008). Evolutionary inference from
Q\textsubscript{ST}.\emph{Molecular Ecology} , \emph{17,} 1885--1896.

Williams, A.P., Cook, E.R., Smerdon, J.E., Cook, B.I., Abatzoglou, J.T.,
Bolles, K., \ldots{} Livneh, B. (2020). Large contribution from
anthropogenic warming to an emerging North American
megadrought.\emph{Science,} \emph{368} , 314-318.

Wright, S. (1931). Evolution in Mendelian populations. \emph{Genetics,
16,} 97--159.

Wright, S.~(1951). The genetic structure of populations.~\emph{Annals of
Eugenics, 15} (1)\emph{,}

323--354.

\textbf{Data Accessibility:}

The trait values from each garden and the SNP dataset used to produce
F\textsubscript{ST} estimates, along with the detailed RADseq protocol
used in this manuscript, will be uploaded to Dryad upon manuscript
acceptance.

\textbf{Author Contributions:}

HF Cooper collected trait data, performed genetic and statistical
analyses, and wrote much of the manuscript. RJ Best contributed to the
development of the Q\textsubscript{ST} plasticity calculation and helped
with the theoretical development of the paper and writing. She, JPM
Corbin, and IJ Garthwaite also assisted with some field collections and
sample processing. LV Andrews helped with RADseq processing. GJ Allan,
CA Gehring, KR Hultine, KR Grady, and TG Whitham, designed the common
garden experiment. All authors provided critical edits to the manuscript
and figures.

\textbf{Tables and Figures:}

\textbf{Table 1:} Model variance and p-values for each trait measured in
each garden, showing the population, genotype, and plot-level effects on
trait variation. In addition, model variance and p-values for each
trait's plasticity across the three gardens, showing the population
effect, are shown.\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
Trait & Garden & Variable & Proportion variance explained &
P-value\tabularnewline
\midrule
\endhead
Bud set & Yuma & Population & 30.19 & \textless{}0.001\tabularnewline
& & Genotype & 13.12 & \textless{}0.001\tabularnewline
& & Plot & 8.43 & \textless{}0.001\tabularnewline
& Agua Fria & Population & 17.70 & \textless{}0.001\tabularnewline
& & Genotype & 24.97 & \textless{}0.001\tabularnewline
& & Plot & 7.82 & \textless{}0.001\tabularnewline
& Canyonlands & Population & 29.53 & \textless{}0.001\tabularnewline
& & Genotype & 16.57 & \textless{}0.001\tabularnewline
& & Plot & 2.23 & 0.003\tabularnewline
& Plasticity & Population & 30.00 & \textless{}0.001\tabularnewline
Bud flush & Yuma & Population & 86.59 & \textless{}0.001\tabularnewline
& & Genotype & 6.42 & \textless{}0.001\tabularnewline
& & Plot & 0 & 1\tabularnewline
& Agua Fria & Population & 78.08 & \textless{}0.001\tabularnewline
& & Genotype & 6.59 & \textless{}0.001\tabularnewline
& & Plot & 1.54 & \textless{}0.001\tabularnewline
& Canyonlands & Population & 2.97 & 0.16\tabularnewline
& & Genotype & 3.56 & 0.08\tabularnewline
& & Plot & 4.07 & 0.02\tabularnewline
& Plasticity & Population & 23.25 & \textless{}0.001\tabularnewline
SLA & Yuma & Population & 24.28 & 0.002\tabularnewline
& & Genotype & 26.81 & \textless{}0.001\tabularnewline
& & Plot & 10.91 & \textless{}0.001\tabularnewline
& Agua Fria & Population & 19.31 & 0.02\tabularnewline
& & Genotype & 8.3 & \textless{}0.001\tabularnewline
& & Plot & 27.63 & \textless{}0.001\tabularnewline
& Canyonlands & Population & 61.21 & \textless{}0.001\tabularnewline
& & Genotype & 8.16 & \textless{}0.001\tabularnewline
& & Plot & 7.03 & \textless{}0.001\tabularnewline
& Plasticity & Population & 6.72 & 0.2\tabularnewline
Height & Yuma & Population & 14.85 & \textless{}0.001\tabularnewline
& & Genotype & 9.40 & \textless{}0.001\tabularnewline
& & Plot & 7.57 & \textless{}0.001\tabularnewline
& Agua Fria & Population & 3.73 & 0.5\tabularnewline
& & Genotype & 11.55 & \textless{}0.001\tabularnewline
& & Plot & 20.26 & \textless{}0.001\tabularnewline
& Canyonlands & Population & 16.39 & 0.004\tabularnewline
& & Genotype & 22.81 & \textless{}0.001\tabularnewline
& & Plot & 8.57 & \textless{}0.001\tabularnewline
& Plasticity & Population & 45.51 & \textless{}0.001\tabularnewline
DRC & Yuma & Population & 12.13 & \textless{}0.001\tabularnewline
& & Genotype & 6.13 & \textless{}0.001\tabularnewline
& & Plot & 12.92 & \textless{}0.001\tabularnewline
& Agua Fria & Population & 0.52 & 0.9\tabularnewline
& & Genotype & 7.90 & \textless{}0.001\tabularnewline
& & Plot & 23.77 & \textless{}0.001\tabularnewline
& Canyonlands & Population & 0.71 & 0.7\tabularnewline
& & Genotype & 24.66 & \textless{}0.001\tabularnewline
& & Plot & 4.03 & \textless{}0.001\tabularnewline
& Plasticity & Population & 36.93 & \textless{}0.001\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2:} Q\textsubscript{ST} and broad-sense heritability,
H\textsuperscript{2} (+ 95\% confidence intervals), for each trait in
each garden, and the mean Q\textsubscript{ST} and H\textsuperscript{2}(+
mean 95\% confidence intervals) for the 100 plasticity permutations for
each trait across gardens.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Trait & Garden & Q\textsubscript{ST} &
H\textsuperscript{2}\tabularnewline
\midrule
\endhead
Bud set & Yuma & 0.54 (0.40-0.86) & 0.19 (0.05-0.26)\tabularnewline
& Agua Fria & 0.26 (0.15-0.42) & 0.30 (0.24-0.34)\tabularnewline
& Canyonlands & 0.47 (0.46-0.57) & 0.28 (0.18-0.24)\tabularnewline
& Plasticity & 0.44 (0.18-0.82) & 0.12 (0.03-0.21)\tabularnewline
Bud flush & Yuma & 0.87 (0.69-0.93) & 0.48 (0.36-0.60)\tabularnewline
& Agua Fria & 0.86 (0.67-0.93) & 0.30 (0.16-0.43)\tabularnewline
& Canyonlands & 0.29 (0.00-1.00) & 0.04 (0.00-0.07)\tabularnewline
& Plasticity & 0.84 (0.64-0.95) & 0.16 (0.04-0.27)\tabularnewline
SLA & Yuma & 0.31 (0.01-0.62) & 0.35 (0.18-0.46)\tabularnewline
& Agua Fria & 0.54 (0.23-0.79) & 0.10 (0.04-0.14)\tabularnewline
& Canyonlands & 0.79 (0.70-0.93) & 0.21 (0.09-0.23)\tabularnewline
& Plasticity & 0.69 (0.37-0.98) & 0.14 (0.01-0.26)\tabularnewline
Height & Yuma & 0.44 (0.23-0.69) & 0.11 (0.04-0.17)\tabularnewline
& Agua Fria & 0.14 (0.00-0.50) & 0.12 (0.05-0.19)\tabularnewline
& Canyonlands & 0.26 (0.03-0.43) & 0.27 (0.18-0.31)\tabularnewline
& Plasticity & 0.66 (0.42-0.94) & 0.09 (0.01-0.17)\tabularnewline
DRC & Yuma & 0.50 (0.09-0.73) & 0.07 (0.04-0.10)\tabularnewline
& Agua Fria & 0.03 (0.00-0.26) & 0.08 (0.05-0.10)\tabularnewline
& Canyonlands & 0.01 (0.00-0.07) & 0.21 (0.19-0.30)\tabularnewline
& Plasticity & 0.07 (0.00-0.16) & 0.18 (0.13-0.23)\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 3:} The adjusted R\textsuperscript{2} and p-value output
from the linear regression models of provenance climate (PC1 score) and
population trait means at each common garden, as well as between PC1 and
the population-level trait plasticity.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
Trait & Garden & Adjusted R\textsuperscript{2} & P-value\tabularnewline
\midrule
\endhead
Bud set & Yuma & 0.67 & \textless{} 0.001\tabularnewline
& Agua Fria & 0.76 & \textless{} 0.001\tabularnewline
& Canyonlands & 0.77 & \textless{} 0.001\tabularnewline
& Plasticity & 0.65 & \textless{} 0.001\tabularnewline
Bud flush & Yuma & 0.49 & 0.002\tabularnewline
& Agua Fria & 0.66 & \textless{} 0.001\tabularnewline
& Canyonlands & -0.01 & 0.36\tabularnewline
& Plasticity & 0.62 & \textless{} 0.001\tabularnewline
SLA & Yuma & -0.08 & 0.69\tabularnewline
& Agua Fria & 0.02 & 0.28\tabularnewline
& Canyonlands & 0.20 & 0.17\tabularnewline
& Plasticity & 0.33 & 0.1\tabularnewline
Height & Yuma & -0.02 & 0.42\tabularnewline
& Agua Fria & 0.32 & 0.013\tabularnewline
& Canyonlands & 0.12 & \textless{} 0.001\tabularnewline
& Plasticity & 0.32 & 0.017\tabularnewline
DRC & Yuma & 0.62 & \textless{} 0.001\tabularnewline
& Agua Fria & -0.01 & 0.37\tabularnewline
& Canyonlands & 0.05 & 0.21\tabularnewline
& Plasticity & 0.28 & 0.02\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 4:} Summary of the two tests for selection. Light blue
indicates a significant Q\textsubscript{ST}-F\textsubscript{ST} test and
dashed lines indicate significant climate regressions for that trait or
plasticity. The Garden abbreviations are Y= Yuma, A = Agua Fria, C =
Canyonlands, and P = Plasticity across gardens.\selectlanguage{english}
\begin{longtable}[]{@{}lllllllllllllllllllll@{}}
\toprule
\textbf{Trait} & Bud set & Bud set & Bud set & Bud set & Bud flush & Bud
flush & Bud flush & Bud flush & SLA & SLA & SLA & SLA & Height & Height
& Height & Height & DRC & DRC & DRC & DRC\tabularnewline
\midrule
\endhead
\textbf{Garden} & Y & A & C & P & Y & A & C & P & Y & A & C & P & Y & A
& C & P & Y & A & C & P\tabularnewline
\textbf{Selection Tests} & & & & & & & & & & & & & & & & & & &
&\tabularnewline
\bottomrule
\end{longtable}

\textbf{Figure 1:} Stand-level mortality event of Fremont cottonwoods
along the Bill Williams National Wildlife Refuge on the lower Colorado
River. Photo taken by HF Cooper in March 2017.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 2:} Map of the 16 collection locations (white circles)
and three common gardens (white stars). The middle common garden of Agua
Fria is also a collection site. Color represents the maximum temperature
of the warmest month (\selectlanguage{ngerman}°C). Axes are labeled with degrees latitude and
longitude. Inset is the map of the United States with a purple rectangle
outlining the latitude and longitude presented in the full map.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 3:} Population means (+/- 1SE) for (a) leaf traits and
(b) growth traits measured in each garden as well as the population mean
trait plasticity values, regressed onto their home climate (indicated by
principal component axis 1 values). Populations are colored by mean
annual temperature (MAT \selectlanguage{ngerman}°C). Regression lines are present when there is
a significant relationship (p \textless{} 0.05) between the PC1 axis and
the trait. SLA is in unit of mm\textsuperscript{2}/mg. Note the scale on
the y-axis for the plasticity regressions are different than the trait
regressions.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 4:} Q\textsubscript{ST} means and 95\% confidence
intervals (point and vertical lines, respectively) for the phenotypic
traits measured at each of the three gardens. The average pairwise
F\textsubscript{ST} value (0.175) + 95\% confidence interval (0.144 -
0.205) is shown as the grey band. Common gardens are abbreviated as Y =
Yuma, A = Agua Fria, and C = Canyonlands.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure 5:} Trait plasticity
Q\textsubscript{ST}-F\textsubscript{ST} distributions for each trait
measured across the three common gardens. The filled shape is the mean
Q\textsubscript{ST} distribution and the open shapes are the 95\%
confidence intervals for the 100 plasticity datasets. Like Fig. 4, the
horizontal lines represent the average population pairwise
F\textsubscript{ST} and its 95\% confidence interval. Any
Q\textsubscript{ST} distribution that crosses the
F\textsubscript{ST}confidence interval is interpreted as no different
from neutral expectation and not under selection.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

