Original Research ARTICLE
Partitioning of multivariate phenotypes using regression trees reveals complex patterns of adaptation to climate across the range of black cottonwood (Populus trichocarpa)
- Department of Forest Resources and Environmental Conservation, Virginia Polytechnic Institute and State University, Blacksburg, VA, USA
Local adaptation to climate in temperate forest trees involves the integration of multiple physiological, morphological, and phenological traits. Latitudinal clines are frequently observed for these traits, but environmental constraints also track longitude and altitude. We combined extensive phenotyping of 12 candidate adaptive traits, multivariate regression trees, quantitative genetics, and a genome-wide panel of SNP markers to better understand the interplay among geography, climate, and adaptation to abiotic factors in Populus trichocarpa. Heritabilities were low to moderate (0.13–0.32) and population differentiation for many traits exceeded the 99th percentile of the genome-wide distribution of FST, suggesting local adaptation. When climate variables were taken as predictors and the 12 traits as response variables in a multivariate regression tree analysis, evapotranspiration (Eref) explained the most variation, with subsequent splits related to mean temperature of the warmest month, frost-free period (FFP), and mean annual precipitation (MAP). These grouping matched relatively well the splits using geographic variables as predictors: the northernmost groups (short FFP and low Eref) had the lowest growth, and lowest cold injury index; the southern British Columbia group (low Eref and intermediate temperatures) had average growth and cold injury index; the group from the coast of California and Oregon (high Eref and FFP) had the highest growth performance and the highest cold injury index; and the southernmost, high-altitude group (with high Eref and low FFP) performed poorly, had high cold injury index, and lower water use efficiency. Taken together, these results suggest variation in both temperature and water availability across the range shape multivariate adaptive traits in poplar.
Tree species and populations may respond to climate change through migration, phenotypic plasticity, evolutionary adaptation, or local extinction (Aitken et al., 2008). Migration rates for trees—estimated to be less than 100 m per year (McLachlan et al., 2005)—are unlikely to be sufficient to keep pace with predicted rates of anthropogenic climate change (IPCC, 2007; Loarie et al., 2009), and adaptation to novel climatic conditions or phenotypic plasticity will, therefore, be necessary for persistence of local populations. Standing genetic variation in traits related to local adaptation increases the likelihood that such populations can make adjustments to mean phenotypes that enables persistence. Indeed, numerous studies have demonstrated that adaptive genetic responses to environmental change can occur in a relatively short period (a few years to hundreds of years), suggesting that local adaptation, at least in part, arises from standing variation rather than new mutations (Savolainen et al., 2013).
Understanding patterns of local adaptation across a species' range facilitates predictions of potential responses to climate change. Forest tree provenance trials that evaluate phenotypic differences among populations have identified traits under divergent selection and quantified their genetic bases (Sork et al., 2013). Such studies allow the (i) exploration of genotype by environment interactions (GxE), which is a necessary condition of local adaptation (Morgenstern, 1996); (ii) comparison of population differentiation in quantitative traits (QST) and neutral molecular markers (FST), which can indicate the role of divergent selection in shaping adaptive genetic differences (Slavov and Zhelev, 2010); and (iii) detection of gradual phenotypic changes along geographical or environmental gradients (clinal variation), which is a signature of local adaptation.
Populus trichocarpa Torr. and Gray (black cottonwood) is a deciduous broadleaf tree species native to western North America that occupies riparian corridors with diverse topography and climate from California to southern Alaska (Gornall and Guy, 2007). Throughout most of its range, P. trichocarpa occurs close to the Pacific coast (within ~200 km), though it extends farther east in the Sierra Nevada Mountains of California and in northern British Columbia. The rapid growth of P. trichocarpa, ease of propagation, intraspecific variation, and substantial hybrid vigor for many important commercial and adaptive traits have contributed to its use for pulp, paper, and bioenergy (Zsuffa et al., 1996), and its potential use for carbon sequestration (Marron et al., 2005). Local adaptation in Populus spp. reflects the integration of many environmental factors—temperature, precipitation, photoperiod, wind, soil nutrient availability, growing season length, and biotic agents—which may be related to geographic variables such as latitude and altitude (Benowicz et al., 2000; Andersson and Fedorkov, 2004; Gornall and Guy, 2007; Hall et al., 2007; Vitasse et al., 2009; Keller et al., 2011; McKown et al., 2013). The strongest evidence of local adaptation in boreal and temperate forest trees is the synchronization of phenotypic traits with local photoperiod and seasonal temperature regimes (Aitken et al., 2008). However, evidence for local adaptation has been found for numerous other traits, including photosynthesis, branch characteristics, nitrogen-use efficiency, and water-use efficiency (Δ13C) (Gornall and Guy, 2007; McKown et al., 2013). Optimization of these traits in a particular environment may involve significant tradeoffs. Understanding these tradeoffs is a fundamental goal of ecological genetics, which may be advanced through multivariate analyses of phenotypes. While clustering and ordination provide homogeneous groupings with respect to multivariate phenotypes, they do not provide direct information on the environmental factors that drive these groupings. Multivariate regression trees (MRT) are an extension of classification and regression trees (CART) that enable such inference through the partitioning of response variables (e.g., phenotypic traits) according to multiple predictors (e.g., climate variables) (Hamann et al., 2011). In this study, we investigated genetic variation in 12 morphological, phenological, and physiological traits, and used MRTs to group these variables according to geography and climate to infer patterns of local adaptation across the range of P. trichocarpa. We also combine estimates of among-population trait differentiation (QST) with a null distribution of genome-wide neutral differentiation (FST) to illustrate the potential role of divergent selection in shaping patterns of trait differentiation.
Materials and Methods
Plant Materials and Common Garden
In 2010, branch cuttings of P. trichocarpa genotypes that span much of the species' range were collected and used to produce plantlets. The plantlets were grown for 6 months in a mist house and subsequently transplanted in May 2011 to an outdoor common garden located at the Reynolds Homestead Forest Resource Research Center located in Patrick County, Virginia (36°37′N and 80°09′W). Ramets of each genotype were planted in a randomized complete block design with 4 blocks. Annual climate variables (means of years 1981–2009) for the closest weather station (Patrick County, VA) showed a mean annual precipitation (MAP) of 1275 mm, a mean annual temperature (MAT) of 14.1°C, a mean warmest month temperature (MWMT) of 24.4°C, and a mean coldest month temperature (MCMT) of 4.2°C (obtained from the weatherbase.com). For this study, we selected 124 genotypes from the 789 available genotypes that covered a wide latitudinal and altitudinal gradient (37–58°N latitude and sea level through ~2300 m) (Figure 1).
Figure 1. Distribution of Populus trichocarpa (green shaded area) and origins of accessions described in this study (red points).
Carbon Isotope Ratio
One branch-section from the current growth year was collected in March 2013 for each tree, oven-dried at 65°C until a constant weight was reached, and ground to 0.5 mm size. About 2 mg of each sample was combusted in an IsoPrime100 isotope ratio mass spectrometer (Isoprime Ltd, Cheadle Hulme, UK) located at the Virginia Tech Forest Soils and Biogeochemistry Laboratory. The carbon isotope ratio, 13C/12C (δ13C), was calculated relative to the international Pee Dee Belemnite standard (Farquhar et al., 1989) as follows:
where Rsa and Rsd are the 13C:12C ratios of the sample and the standard, respectively (Craig, 1957).
Growth and Branching
In March 2013, while the trees were in a dormant state, growth and branch characteristics (Figure 2) were measured on the 124 genotypes. Height (H), crown diameter (CD), insertion height of the highest branch (proleptic or sylleptic) (Ih), and insertion height of the lowest branch (Il) (Broeckx et al., 2012), were measured using a tape and a telescopic pole. Stem diameter (D) at 22 cm above ground level was measured using a digital caliper. Volume Index (VI = D2xH) was calculated according to Causton (1985). We also counted the total number of branches (NB) and the number of sylleptic branches produced during the 2012 growing season (NSyll) on the main stem of each tree. The relative number of branches (RNB = density of branches per unit of stem height) was calculated by dividing NB by H. CD was obtained by computing the mean of the crown diameters taken from two perpendicular directions (east to west and north to south). The relative canopy depth (RCD), defined as the percentage of the stem carrying branches, was computed as (Ih–Il)/H.
Figure 2. Illustration of growth and branching traits. H, Height; D, Stem diameter at 22 cm above ground level; CD, Crown Diameter (measured North-South and East-West); Ihigh, insertion height of the highest branch; Ilow, height of the lowest branch.
Bud Phenology and Cold Injury
Bud stage was scored on a weekly basis during April 2012 (BF = bud flush) and September 2012 (BS = bud set), and the date of BF or BS was converted to the Julian day (the number of days from December 31). The bud was considered flushed when a leaf had emerged 1 cm from the bud, and set when dark red/brown scales covered the apex. Buds that set or flushed before the first week or after the last week of scoring were recorded as having occurred 1 week before or 1 week later, respectively. Fall cold injury (I−20) was assessed by measuring electrolytic leakage of branch sections collected on November 5, 2012 (Hannerz et al., 1999). Cool (~4°C) temperatures were recorded at the site prior to sampling for cold hardiness, but no freeze event had yet occurred. For each genotype, three shoot segments (1–2 mm) were placed in solutions of 500 μl of distillated water with a trace amount of silver iodide to facilitate ice nucleation. Control samples were kept at 4°C. Test samples were placed in a Tenney programmable freezing chamber and the temperature was gradually lowered at a rate of 4°C/h until they reached –20°C, at which point they were incubated for 2 h and subsequently thawed at 4°C. An additional 500 μl of distilled water was then added and the samples were agitated on an orbital shaker to homogenize the solution. Conductivity was measured with a digital conductivity meter following the temperature treatments, and after heat-kill (95°C for 4 h). An index of freezing injury (It) was calculated according to Flint et al. (1967):
Where It is the index of injury (or percent injury), Rt and R0 are the relative conductivities for –20°C and 4°C respectively, Lt is the specific conductivity of leachate from the sample frozen at –20°C, Lk is the specific conductivity of leachate from the sample frozen at –20°C and then heat-killed, L0 is the specific conductivity of leachate from control samples (4°C), and Ld is the specific conductivity of leachate from the heat-killed control samples (Flint et al., 1967).
Populations Defined by Multivariate Regression Trees (MRT)
Data were visually inspected for normality and homoscedasticity, and NSyll was linearized by log transformation before analysis. To reduce the number of variables, and to deal with the high correlation between growth and branching variables, principal component analysis (PCA) was run on the growth and branch traits using the prcomp function in R (R Core Team, 2014). MRTs were constructed using the R package MVpart v1.2–6 (http://cran.rproject.org/web/packages/mvpart/index.html). MRT is an extension of Classification and Regression Trees (CART) that can be used to explore, describe, and predict relationships between multivariate response data and multiple predictor variables (De'ath, 2002). A set of clusters are grown by repeated binary splits of the dataset to produce nodes as homogeneous as possible with respect to the response variables. The geographic and climatic variables were used as predictors and two principal components, BS, BF, I–20, and δ13C were used as response variables. For this analysis, 21 annual and seasonal climate variables (means of 1981—2009) were obtained using the ClimateWNA program (Wang et al., 2006) (v4.72). In addition, a drought index, the Standard Precipitation Evapotranspiration Index (SPEI), was obtained from the Global SPEI database (SPEIbase v.2.2). SPEI is calculated from the Climate Research Unit's TS3.2 database (http://www.cru.uea.ac.uk/data). Descriptions of each variable can be found in Supplemental Table 1.
Estimation of Genetic Parameters
Samples were grouped into 20 populations according to their geographic origins (latitude, longitude, and altitude) for computation of population differentiation in quantitative traits (QST) and broad-sense heritability (H2). For each trait, the clonal best linear unbiased predictors (BLUPs) and variance components were estimated using restricted maximum likelihood (REML) with block as a fixed factor and genotype and population as random factors according to the following model:
where Yijk is the value for the kth genotype from the jth population in the ith block, μ is the population mean, Bi is the fixed effect of the ith block, Pj is the random effect of the jth population, Gk(Pj) is the random effect of the kth genotype in the jth population, and εijk is the experimental error. QST was computed as:
Where σ2P is the among population component of variance and σ2G(P) is the genotype-within-population component of variance. These variance component estimates are based on clonal replication of trees, and have been shown by Lynch and Walsh (1998) to give QST values confounded by non-additive genetic variance and maternal effects, which result in a slight deflation of QST. Hall et al. (2007) showed that QST values estimated using this method probably represent lower bounds of QST. Therefore, our values are probably conservative with respect to adaptive population differentiation. Broad-sense heritabilities were calculated across-population (H2G) and within-population (H2G(P)):
Finally, genetic correlations between traits were obtained by calculating correlations among the respective BLUPs, and environmental correlations were estimated according to Searle (1961):
Where R, r, and r′ are respectively the phenotypic, genetic, and environmental correlations, and H2G1 and H2G2 are the across-population heritabilities (H2G) of the relevant traits. Parametric bootstrapping was used to compute the 95% confidence intervals for the genetic parameters (Spitze, 1993; O'Hara and Merila, 2005) with 1000 bootstrap iterations of the variances components.
Genotyping and Population Structure
We previously reported using sequence capture to retrieve and sequence much of the exome for 48 P. trichocarpa individuals originating from across the natural range (Zhou and Holliday, 2012; Zhou et al., 2014). These data have now been augmented with sequencing of an additional 391 samples. Capture, sequencing, assembly, and alignment are detailed in Zhou and Holliday (2012). Briefly, young leaf tissue was collected and genomic DNA was extracted using a Qiagen DNeasy Plant Mini Kit (Qiagen, Inc, Valencia, California, USA). A custom RNA bait library was designed and synthesized by Agilent Technologies (Santa Clara, California, USA), which targeted at least one exon for each gene model, 120 bp upstream of each gene model, as well as 1,000 random intergenic regions (at least 1000 bp from any annotated gene model). Our previous reports used a library comprised of 173,040 baits. Due to changes in the Agilent product, we were able to increase this to 230,720, a ~33% increase in the number of targeted regions. The data reported here includes only samples arising from this revised bait library. Captured libraries were sequenced on an Illumina HiSeq instrument. Following sequencing, the raw reads were trimmed of low quality bases, adapter sequences, demultiplexed using custom scripts, and aligned to the P. trichocarpa reference genome using BWA software (Li and Durbin, 2010). These reads were trimmed to retain the maximum length with an average quality score of 30 and a minimum of 13, with no undetermined nucleotides. For each genotype, SNPs were called for the uniquely mapped and re-paired reads using SAMtools software (Li et al., 2009) using a combined 8x cutoff of coverage depth and quality score of 30. A total of 1,058,513 SNPs were retained from the genomic regions that have 8x coverage in at least 75% of the 391 clones. Using the genotypes for which phenotyping is described in this paper, we estimated mean pairwise FST for the 727,993 segregating SNPs in that sample using the hierfstat R package according to the geographic groupings described above. To identify quantitative traits that may be under divergent selection, we compared QST estimates with the empirical distribution of FST, which was visualized using the smooth.spline function in R with a smoothing parameter of 0.6.
Relationships among Traits, Geography, and Climate
Correlations were computed for each pair of variables to evaluate the relationships among traits, geography, and climate. Among geographical variables, latitude had the highest number of significant correlations with climate variables (20/22 significant and 10/22 with |r| > 0.6) (Table 1). Longitude and altitude also had a moderate correlation with temperature and precipitation variables, but due to the southeast—northwest orientation of the species range, as well as the high altitude of the southernmost samples, it is difficult to disentangle latitude from longitude and altitude. The exception to this was the central part of the range (between 48° and 54°N Lat.), for which there was an increase in temperature differential (TD) and water deficit from the coast to the interior.
Table 1. Pearson correlation coefficients and levels of significance between climate and geographic variables.
All traits except RNB had a significant (p < 0.05) relationship with latitude of origin (Table 2). BS had the highest correlation with latitude (r = –0.79), while δ13C had correlations of 0.20 and –0.46 with latitude and altitude, respectively. The relationship with altitude was mostly driven by the southernmost group of samples, which had a low mean δ13C (–24.77) compared with the other populations. Among the traits, the highest correlations with climate were with temperature related variables (MAT, MWMT, MCMT, degree-days (DD), frost-free period (FF), EMT, EXT; Supplemental Table 2). BS showed the highest correlations with the climate variables, while δ13C had a weak negative but significant correlations with moisture variables (Hargreaves climatic moisture deficit (CMD): r = –0.32; Annual heat:moisture index (AHM): r = –0.22; Summer heat:moisture index (SHM): r = –0.27; Hargreaves reference evaporation (Eref): r = –0.28). Relatively strong positive genetic correlations (r > 0.50) were found between growth traits (H, D, VI) and branch traits (NB, NSyll and CD) (Table 3), with large trees (H, D, VI) generally displaying the highest NB and NSyll. I–20 was correlated with BS (r = 0.53) but there was no correlation between BF and either I–20 or BS. However, growth and branching traits had significant negative relationships with BF, and a positive relationship with BS and I–20. Finally, δ13C had had no significant correlations with the other traits.
Table 2. Pearson correlation coefficients (r) and levels of significance calculated from genotypic means between traits and geographic variables.
Table 3. Genetic (below diagonal) and environmental (above diagonal) correlations between traits with levels of significance.
Multivariate Regression Trees
Prior to MRT analysis, growth and branching traits were subjected to PCA to remove their colinearity. The first two principal components were retained because they explained most of the variation in growth and branching traits (88%). The loadings of the traits (Supplemental Table 3) indicate that the first principal component (PC1) was positively correlated with all eight variables (0.56–0.95) whereas PC2 was mainly reflects branch characteristics (NB, NSyll, RNB, and RCD). The MRT using geography predictors split the 124 genotypes into five groups, which explained 42% of the variation in the traits (Figure 3). Trees originating from below 49.28°N Lat. were split into two groups according to their altitude of origin, a split that could also have been made based on latitude (40.30°N Lat.) or longitude (120.9°W Long.) with the same amount of variance explained. The seven southernmost genotypes (Group 1) had lower than average growth, fewer branches, and the lowest δ13C. This group also flushed early (small BF), set buds late (high BS), and had a high I–20. The 32 genotypes from Group 2, originating between 49.34° and 40.30°N Lat. (coastal), had above average growth and branching, small BF, high BS, and the highest I–20. Trees originating between latitudes 51.05° and 49.28°N Lat. (southern BC) were further differentiated into Groups 3 and 4 according to their longitude of origin. Group 4 (interior group) originated from areas east of 123°W and consisted of 29 genotypes characterized by growth (H, D, CD) slightly below the mean and low I–20, whereas the 24 genotypes from Group 3 (coastal group) originated from areas west of 123°W Long. and had growth and I–20 slightly above the sample mean. Group 5 consisted of 31 genotypes originating above 51.05°N. This group had the lowest growth, fewest branches, latest BF, earliest BS, and lowest I–20.
Figure 3. Multivariate regression tree (MRT) analyses of the 124 black cottonwood genotypes with latitude, longitude and altitude as predictor variables (A); with climate as predictor variables (B); and geographic representation of the splits (C,D). Bar charts at the terminal leaves of the regression trees represent population means for the traits.
MRTs using climate variables as predictors also split the 124 genotypes in five groups, which explained 39% of variation in the measured traits (Figure 3). Group 1 contained seven high altitude genotypes, which originated from areas with Eref greater than 744 and FFP less than 204.5 days. These trees had low growth and branching, early BF, late BS and the lowest δ13C. Groups 2 and 3 originated from areas with Eref greater than 744 and FFP higher than 204.5 days, and were further split according to MAP, which revealed spatial heterogeneity in precipitation for this region. The 19 genotypes from Group 2 (MAP higher than 1004 mm) exhibited early BF, had the latest BS, the highest I–20, the highest performance in growth and branch characteristics, and the highest δ13C. The eight genotypes of Group 3 (MAP < 1004 mm) had growth and branch traits close to the average, the earliest BF, late BS, and high I–20. Group 4 had 29 genotypes originating from an area characterized by Eref < 744 and MWMT > 16.85°C, respectively. Genotypes from this population were characterized by growth and branching traits, δ13C, BS, BF, and I–20 close to the average. Finally, Group 5 was comprised of 31 genotypes that originate from localities with Eref lower than 744 and MWMT higher than 16.85°C. Genotypes from this group had the lowest growth and branch traits, early BS, late BF, and the lowest I–20.
Population Differentiation and Heritability
BS and I–20 had the highest among-population differentiation, with QST ≥ 0.5. Growth traits (H, D and VI) and δ13C had moderate values (0.32 ≤ QST ≤ 0.44), while differentiation for branch traits (NB, NSyll, RNB and RCD) and BF was low (QST < 0.25) (Table 4). BS, I–20, growth traits, and δ13C exceeded the 99th percentile of the empirical distribution of FST, while the branching traits and BF fell below this threshold (Figure 4). Heritabilities (H2G(P)) were low to moderate (0.13 ≤ H2 ≤ 0.32; Table 4). BF, I–20, NB, and NSyll had the highest values (H2G(P) ≥ 0.29), while VI had the lowest (H2 = 0.13).
Table 4. Population differentiation in quantitative traits (QST) and broad sense heritability within-population (H2G(P)) and across the collection (H2G) for the 12 traits.
Figure 4. Empirical distribution of FST based on 727,993 SNPs. Dotted lines indicate the 99th percentile of FST values. QST estimates for each trait are indicated with colored triangles.
Local Adaptation vs. Drift
In spite of high gene flow among populations, temperate and boreal trees often display substantial adaptive genetic differentiation due to strong climatic gradients across their ranges (Farmer, 1996; Morgenstern, 1996; Howe et al., 2003; Savolainen et al., 2007; Aitken et al., 2008). In general, long-lived and widely distributed woody species harbor low levels of population differentiation in molecular markers, with mean FST often less than 0.10 (Hamrick et al., 1992). For Populus, in particular, FSTis even lower (<0.05; Slavov and Zhelev, 2010). Comparisons between FST and QST are problematic when the distribution of the former, which may be influenced by demographic history, is not known (Whitlock, 2008). We, therefore, used a large, genome-wide sample of SNPs to estimate the empirical distribution of FST and determine the relative influence of drift and selection in generating observed patterns of trait differentiation across the range. Each of H, D, VI, CD, δ13C, BS, and I–20 fell outside the 99th percentile of the empirical distribution of FST, which suggests that their differentiation is shaped by divergent natural selection.
Populus species exhibit strong latitudinal clines in bud set (Pauley and Perry, 1954), and the high values of QST we found for bud set and cold injury are in agreement with estimates in range-wide studies of other woody taxa (Morgenstern, 1996; Hurme et al., 1997; Hall et al., 2007; Mimura and Aitken, 2007; Keller et al., 2011). Autumn phenology and cold acclimation are important targets of natural selection in northern tree species because they enable synchronization of annual growth cycle with the local climatic cycle to avoid cold damage while maximizing growth (Howe et al., 2003). Morgenstern (1969) and Keller et al. (2011) found high QST values for BF (QST = 0.52 in Picea mariana and 0.66 in Populus balsamifera), but most studies in woody species have reported low QST for BF, in agreement with our results (Howe et al., 2003). This may reflect low genetic variation for the trait, but also the nature of many common garden studies, in which populations from across a latitudinal gradient are grown in an intermediate location. As bud flush is governed by both a chilling and heat sum requirement, it may be that when grown in such intermediate locations, trees from the north meet their heat-sum requirements earlier than they do at their origin (spring temperatures warmer than the population origin), whereas trees from the south meet that requirement later than at their origin (spring temperatures colder than the population origin). This hypothesis assumes that all populations receive adequate chilling. The high QST value for bud set in P. balsamifera (Keller et al., 2011), the sister species to P. trichocarpa, may reflect the location of the common garden in that study, which was at the southern limit of the species range and included samples from across ~15° of latitude. Where among-population genetic differentiation in chilling or heat-sum requirements exists, such a common garden location may better reveal these differences. Further research using reciprocal transplants is needed to test this hypothesis. Finally, while height and diameter also exhibited strong differentiation consistent with local adaptation, QST for branch traits fell below our empirical FST threshold, suggesting that any divergent selection related to branch traits is weak.
Patterns of Adaptation to Local Climate
Genotypes from the northern periphery of our sampling area had the lowest growth, which probably reflects adaptation to a short growing season (shorter FFP) in their native environment. For many tree species, common garden experiments have revealed that high latitude populations achieve less height growth even when they display higher assimilation rates than low latitude populations (Burtt, 1956; Kaurin et al., 1985; Junttila and Kaurin, 1990; Soolanayakanahally et al., 2009; Savage and Cavender-Bares, 2013). Interestingly, the southernmost genotypes (originating from between 1227 and 2363 m above sea level) had similar growth to those from the north, which could be explained by adaptation to the high altitude climate conditions in the Sierra Nevada Mountains of California. This area has low minimum temperatures and a short growing season (median EMT = –22.5 and FFP = 169.3) similar to the north (EMT = –25.8 and FFP = 179), and is much cooler than the common garden site. Climate variables related to water availability also reveal drier summers for the Sierra Nevada location (MSP = 159.2, Eref = 601, and CMD = 100) compared to the rest of the species range (mean rangewide MSP = 290.1, Eref = 728.8 and CMD = 273.3). Decreasing tree height with increasing altitude and latitude of origin has been described for a number of conifers and broadleaf trees (Rehfeldt, 1994; Oleksyn et al., 1998; Sáenz-Romero et al., 2006; Premoli et al., 2007; Rweyongeza et al., 2007; Vitasse et al., 2009). In addition to the tradeoff between growth and cold hardiness, possible explanations for the low growth of high altitude populations include high respiration rates, high allocation to roots, or reduced shoot-growth period (Oleksyn et al., 1998). At the same time, the high altitude genotypes from California had a relatively long growing season (early BF and late BS) in the common garden in spite of having a lower FFP (169 days) in their locality of origin. These two environments differ in temperature regime but not significantly in photoperiod (only 1° latitude), which suggests that the long growth period that occurred in the common garden was likely related to temperature differences. Indeed, the higher temperatures in the common garden compared to high altitude locations may have led to an earlier BF, which is mainly controlled by temperature (Luquez et al., 2008; Keller et al., 2011), although insufficient chilling may delay BF.
The center of the range (southern BC, Oregon, Washington) was partitioned into three populations. Genotypes from southern BC, which is characterized by moderate temperatures, high precipitation, and low evaporative demand, were differentiated into two groups according to continentality: the coastal group, which originates from a mild climate (MWMT > 16.85) with a long growing season, displayed better growth (slightly above the overall mean), and was more susceptible to fall cold damage than the interior group. The interior population in southern BC was merged with the northern groups in the MRT using climate predictors, as the temperature variables of these two population locations were similar. Genotypes spanning the coast of Oregon and Washington were the tallest, and had the highest number of branches in the common garden. This area is below 400 m altitude and has a mild climate (highest MAT, MWMT, EXT and low TD) that is similar to the common garden in terms of temperature variables. The climate variable-based MRT showed that these genotypes could be further differentiated according to MAP in spite of closely overlapping geography. This local spatial heterogeneity in climate could help these populations respond to climate change, facilitating migration by reducing dispersal distance required to find an appropriate climatic niche (Ackerly et al., 2010).
Plants adapted to arid environments generally have higher δ13C, an indirect proxy for water use efficiency (WUE), than plants adapted to wet environments (Passioura, 1982). By contrast, we found the lowest δ13C for genotypes originating from the most water-limited areas. Similar trends were observed in Pseudotsuga menziesii, Larix occidentalis, Pinus ponderosa, and Pinus contorta (Zhang et al., 1993; Zhang and Marshall, 1994, 1995; Aitken et al., 1995; Guy and Holowachuk, 2001), and may reflect differences in environmental conditions between common garden and population origins (Read and Farquhar, 1991; Aitken et al., 1995). Plants cope with fluctuation in water availability by adjusting their leaf area, root depth and density, hydraulic properties, photosynthetic capacity, and stomatal conductance. The adaptive strategy they use determines the strength and the sign of the correlation between δ13C and the climate (Aitken et al., 1995; Lauteri et al., 2004). The Sierra Nevada Mountains population (Group 1), which had the lowest δ13C, originates from a climate characterized by drier summers (MSP = 159 mm, AHM = 19.0, SHM = 135.8, Eref = 991) compared with the other populations (mean MSP = 328 mm, AHM = 13.3, SHM = 64.2, Eref = 681.5), for which δ13C values were higher. This may seem counterintuitive as one may expect genotypes adapted to dry locations would exhibit higher WUE (i.e., higher δ13C). However, when grown in the common garden, genotypes of Group 1 likely encountered a moister summer than in their native environment, and consequently may have adjusted one of the physiological parameters noted above, which may have resulted in lower δ13C. On the other hand, accessions from populations originating from areas with greater summer precipitation may have experienced mild water stress and responded by reducing their water loss, which led to higher δ13C. Aitken et al. (1995) found a similar result for douglas-fir (Pseudotsuga menziesii) seedlings grown at a location with intermediate precipitation compared with source populations (coastal and interior). This hypothesis is also supported by δ13C lower than the average for Group 3 in the climate-based MRT, for which summer moisture values (MSP = 138 mm, AHM = 21.4, SHM = 120.5, Eref = 865) were similar to Group 1.
Relationships Between Productivity and Adaptability Traits
Branch characteristics are closely related to productivity as they determine the quantity of light interception and CO2 assimilation (Halle et al., 1978), which may explain the positive correlation between growth (H, D, VI) and branching variables. We found a higher number of sylleptic branches on highly productive genotypes, which is consistent with Ceulemans et al. (1990) and Scarascia-Mugnozza et al. (1999), who showed that sylleptic branches supply a larger proportion of carbon to the stem than proleptic branches in young poplar trees. The lack of a negative correlation between δ13C and growth/branch characteristics combined with moderate H2 for these traits suggests that, for P. trichocarpa, there is a possibility of improving WUE without compromising productivity. The continued development of poplars as bioenergy feedstocks depends partially on their suitability for cultivation on marginal lands that may be water limited (Berndes, 2002). Such a strategy would limit competition with food crops for the best sites. Hence, WUE will be a key target of selection of improved cultivars. Similar results were reported for interspecific poplar hybrids (Rae et al., 2004; Marron et al., 2005; Monclus et al., 2005, 2006; Fichot et al., 2010), P. nigra (Chamaillard et al., 2011), and P. trichocarpa (McKown et al., 2013). Farquhar et al. (1989) suggested that a lack of relationship between δ13C and growth means that stomatal conductance drives WUE variation, while a positive correlation between the two traits generally supports that differences in photosynthetic capacity is the main driver for WUE. Gornall and Guy (2007) showed that differences among P. trichocarpa genotypes in WUE were related to variability in stomatal conductance rather than assimilation rate.
Our data suggest that P. trichocarpa has considerable standing adaptive variation within and among populations, which represent an adaptive potential. We found strong population differentiation in growth, phenology, cold hardiness, and carbon isotope composition, but weak genome-wide differentiation at SNP loci, which suggests that local adaptation better explains patterns of variation in these traits than drift alone. Multivariate regression trees using geography or climate variables as predictors showed that population differentiation in adaptive traits could be primarily explained by either latitude or evapotranspiration, respectively. The MRT analysis also suggested that the geographically intermingled accessions from coastal California and Oregon could be separated for the traits we measured on the basis of heterogeneity in mean annual precipitation. Genotypes from these areas were the best adapted to the climatic conditions of our Virginia common garden, suggesting that they can be used for breeding purposes in areas with similar environmental conditions. Indeed, with the growing interest in poplar species for bioenergy, knowledge of suitable genotypes for specific environments constitutes valuable information for successful plantations. One limitation of this study is the location of the common garden, which was located well outside the natural longitudinal range of P. trichocarpa, and at a latitude similar to the southern species range limit. In addition to the impacts of daylength regimes the common garden location likely had on growth and phenology, this environment may also impact the rankings of genotypes due to differences in soil conditions and temperature. In addition, fungal pathogens present in the southeastern US may impact performance of this species, particularly Septoria musiva, though we did not observe cankers during the course of this study. To address these issues, we have established an additional common garden site near the center of the species range, which will help to tease apart the relationship between genetics and environment.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We would like to thank Kyle Peer, Clay Sawyers, and Deborah Bird (Virginia Tech Reynold's Homestead Forestry Research Station) for assistance with plant propagation, establishment, and maintenance of the common garden, as well as the following organizations and individuals for providing genetic material for this study: Drs. Chang-Yi Xie and Alvin Yanchuk (British Columbia Ministry of Forests, Lands, and Natural Resource Operations), Dr. Brian Stanton (Greenwood Resources), Dr. Brad St. Clair [United States Forest Service (USFS)], and Dr. Dennis Ringes (USFS, retired). Finally, we thank two reviewers whose insightful comments greatly improved this manuscript. This work was supported by the National Science Foundation Plant Genome Research Program (IOS: 1054444) (grant to JH) and a Fulbright Institute of International Education Fellowship to RWO.
The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fpls.2015.00181
Ackerly, D., Loarie, S., Cornwell, W., Weiss, S., Hamilton, H., Branciforte, R., et al. (2010). The geography of climate change: implications for conservation biogeography. Divers. Distrib. 16, 476–487. doi: 10.1111/j.1472-4642.2010.00654.x
Aitken, S. N., Kavanagh, K. L., and Yoder, B. J. (1995). Genetic variation in seedling water-use efficiency as estimated by carbon isotope ratios and its relationship to sapling growth in Douglas-fir. For. Genet. 2, 199–206.
Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T., and Curtis-Mclane, S. (2008). Adaptation, migration or extirpation: climate change outcomes for tree populations. Evol. Appl. 1, 95–111. doi: 10.1111/j.1752-4571.2007.00013.x
Benowicz, A., Guy, R. D., and El-Kassaby, Y. A. (2000). Geographic pattern of genetic variation in photosynthetic capacity and growth in two hardwood species from British Columbia. Oecologia 123, 168–174. doi: 10.1007/s004420051002
Broeckx, L. S., Verlinden, M. S., Vangronsveld, J., and Ceulemans, R. (2012). Importance of crown architecture for leaf area index of different Populus genotypes in a high-density plantation. Tree Physiol. 32, 1214–1226. doi: 10.1093/treephys/tps083
Causton, D. (1985). “Biometrical, structural and physiological relationships among tree parts,” in Attributes of Trees as Crop Plants, eds M. G. R. Cannel and J. E. Jackson (Cumbria: Titus Wilson and Son Ltd.), 137–159.
Ceulemans, R., Stettler, R. F., Hinckley, T. M., Isebrands, J. G., and Heilman, P. E. (1990). Crown architecture of Populus clones as determined by branch orientation and branch characteristics. Tree Physiol. 7, 157–167. doi: 10.1093/treephys/7.1-2-3-4.157
Chamaillard, S., Fichot, R., Vincent-Barbaroux, C., Bastien, C., Depierreux, C., Dreyer, E., et al. (2011). Variations in bulk leaf carbon isotope discrimination, growth and related leaf traits among three Populus nigra L. populations. Tree Physiol. 31, 1076–1087. doi: 10.1093/treephys/tpr089
Craig, H. (1957). Isotopic standards for carbon and oxygen and correction factors for mass-spectrometric analysis of carbon dioxide. Geochim. Cosmochim. Acta 12, 133–149. doi: 10.1016/0016-7037(57)90024-8
Farmer, R. E. Jr. (1996). “The genecology of Populus,” in Biology of Populus and its Implications for Management and Conservation, eds R. F. Stettler, H. D. Bradshaw, P. E. Heilman, and T. M. Hinckley (Ottawa, ON: NRC Research Press), 33–55.
Farquhar, G. D., Ehleringer, J. R., and Hubick, K. T. (1989). Carbon isotope discrimination and photosynthesis. Annu. Rev. Plant Physiol. Plant Mol. Biol. 40, 503–537. doi: 10.1146/annurev.pp.40.060189.002443
Fichot, R., Barigah, T. S., Chamaillard, S., Thiec D. L. E., Laurans, F., Cochard, H., et al. (2010). Common trade-offs between xylem resistance to cavitation and other physiological traits do not hold among unrelated Populus deltoides x Populus nigra hybrids. Plant Cell Environ. 33, 1553–1568. doi: 10.1111/j.1365-3040.2010.02164.x
Flint, H. L., Boyce, B. R., and Beattie, D. J. (1967). Index of injury—a useful expression of freezing injury to plant tissues as determined by the electrolytic method. Can. J. Plant Sci. 47, 229–230. doi: 10.4141/cjps67-043
Guy, R. D., and Holowachuk, D. L. (2001). Population differences in stable carbon isotope ratio of Pinus contorta Dougl. ex Loud.: relationship to environment, climate of origin, and growth potential. Can. J. Bot. 79, 274–283. doi: 10.1139/b01-001
Hall, D., Luquez, V., Garcia, V. M., St Onge, K. R., Jansson, S., and Ingvarsson, P. K. (2007). Adaptive population differentiation in phenology across a latitudinal gradient in European aspen (Populus tremula, L.): a comparison of neutral markers, candidate genes and phenotypic traits. Evolution 61, 2849–2860. doi: 10.1111/j.1558-5646.2007.00230.x
Howe, G. T., Aitken, S. N., Neale, D. B., Jermstad, K. D., Wheeler, N. C., and Chen, T. H. (2003). From genotype to phenotype: unraveling the complexities of cold adaptation in forest trees. Can. J. Bot. 81, 1247–1266. doi: 10.1139/b03-141
IPCC. (2007). Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge; New York, NY: Cambridge University Press.
Keller, S. R., Soolanayakanahally, R. Y., Guy, R. D., Silim, S. N., Olson, M. S., and Tiffin, P. (2011). Climate-driven local adaptation of ecophysiology and phenology in balsam poplar, Populus balsamifera L. (Salicaceae). Am. J. Bot. 98, 99–108. doi: 10.3732/ajb.1000317
Lauteri, M., Pliura, A., Monteverdi, M. C., Brugnoli, E., Villani, F., and Eriksson, G. (2004). Genetic variation in carbon isotope discrimination in six European populations of Castanea sativa Mill. originating from contrasting localities. J. Evol. Biol. 17, 1286–1296. doi: 10.1111/j.1420-9101.2004.00765.x
Luquez, V., Hall, D., Albrectsen, B., Karlsson, J., Ingvarsson, P., and Jansson, S. (2008). Natural phenological variation in aspen (Populus tremula): the SwAsp collection. Tree Genet. Genomes 4, 279–292. doi: 10.1007/s11295-007-0108-y
Marron, N., Villar, M., Dreyer, E., Delay, D., Boudouresque, E., Petit, J. M., et al. (2005). Diversity of leaf traits related to productivity in 31 Populus deltoides x Populus nigra clones. Tree Physiol. 25, 425–435. doi: 10.1093/treephys/25.4.425
McKown, A. D., Guy, R. D., Klapste, J., Geraldes, A., Friedmann, M., Cronk, Q. C., et al. (2013). Geographical and environmental gradients shape phenotypic trait variation and genetic structure in Populus trichocarpa. New Phytol. 201, 1263–1276. doi: 10.1111/nph.12601
Monclus, R., Dreyer, E., Delmotte, F. M., Villar, M., Delay, D., Boudouresque, E., et al. (2005). Productivity, leaf traits and carbon isotope discrimination in 29 Populus deltoides x P. nigra clones. New Phytol. 167, 53–62. doi: 10.1111/j.1469-8137.2005.01407.x
Monclus, R., Dreyer, E., Villar, M., Delmotte, F. M., Delay, D., Petit, J.-M., et al. (2006). Impact of drought on productivity and water use efficiency in 29 genotypes of Populus deltoides x Populus nigra. New Phytol. 169, 765–777. doi: 10.1111/j.1469-8137.2005.01630.x
Oleksyn, J., Modrzýnski, J., Tjoelker, M. G., Z·Ytkowiak, R., Reich, P. B., and Karolewski, P. (1998). Growth and physiology of Picea abies populations from elevational transects: common garden evidence for altitudinal ecotypes and cold adaptation. Funct. Ecol. 12, 573–590. doi: 10.1046/j.1365-2435.1998.00236.x
Passioura, J. B. (1982). “Water in the soil-plant-atmosphere continuum,” in Physiological Plant Ecology II, eds O. L. Lange, P. S. Nobel, C. B. Osmond, and H. Ziegler (Berlin; Heidelberg: Springer), 5–33.
Premoli, A. C., Raffaele, E., and Mathiasen, P. (2007). Morphological and phenological differences in Nothofagus pumilio from contrasting elevations: evidence from a common garden. Austral Ecol. 32, 515–523. doi: 10.1111/j.1442-9993.2007.01720.x
Rae, A. M., Robinson, K. M., Street, N. R., and Taylor, G. (2004). Morphological and physiological traits influencing biomass productivity in short-rotation coppice poplar. Can. J. For. Res. 34, 1488–1498. doi: 10.1139/x04-033
Sáenz-Romero, C., Guzmán-Reyna, R. R., and Rehfeldt, G. E. (2006). Altitudinal genetic variation among Pinus oocarpa populations in Michoacán, Mexico: implications for seed zoning, conservation, tree breeding and global warming. For. Ecol. Manage. 229, 340–350. doi: 10.1016/j.foreco.2006.04.014
Savage, J. A., and Cavender-Bares, J. (2013). Phenological cues drive an apparent trade-off between freezing tolerance and growth in the family Salicaceae. Ecology 94, 1708–1717. doi: 10.1890/12-1779.1
Scarascia-Mugnozza, G., Hinckley, T., Stettler, R., Heilman, P., and Isebrands, J. (1999). Production physiology and morphology of Populus species and their hybrids grown under short rotation. III. Seasonal carbon allocation patterns from branches. Can. J. For. Res. 29, 1419–1432. doi: 10.1139/cjfr-29-9-1419
Slavov, G. T., and Zhelev, P. (2010). “Salient biological features, systematics, and genetic variation of Populus,” in Genetics and Genomics of Populus, eds S. Jansson, R. Bhalerao, and A. Groover (New York, NY: Springer), 15–38. doi: 10.1007/978-1-4419-1541-2_2
Soolanayakanahally, R. Y., Guy, R. D., Silim, S. N., Drewes, E. C., and Schroeder, W. R. (2009). Enhanced assimilation rate and water use efficiency with latitude through increased photosynthetic capacity and internal conductance in balsam poplar (Populus balsamifera L.). Plant Cell Environ. 32, 1821–1832. doi: 10.1111/j.1365-3040.2009.02042.x
Sork, V., Aitken, S., Dyer, R., Eckert, A., Legendre, P., and Neale, D. (2013). Putting the landscape into the genomics of trees: approaches for understanding local adaptation and population responses to changing climate. Tree Genet. Genomes, 9, 901–911. doi: 10.1007/s11295-013-0596-x
Vitasse, Y., Delzon, S., Bresson, C. C., Michalet, R., and Kremer, A. (2009). Altitudinal differentiation in growth and phenology among populations of temperate-zone tree species growing in a common garden. Can. J. For. Res. 39, 1259–1269. doi: 10.1139/X09-054
Wang, T., Hamann, A., Spittlehouse, D. L., and Aitken, S. N. (2006). Development of scale-free climate data for western Canada for use in resource management. Int. J. Climatol. 26, 383–397. doi: 10.1002/joc.1247
Zhang, J., and Marshall, J. (1995). Variation in carbon isotope discrimination and photosynthetic gas exchange among populations of Pseudotsuga menziesii and Pinus ponderosa in different environments. Funct. Ecol. 9, 402–412. doi: 10.2307/2390003
Zhou, L., Bawa, R., and Holliday, J. A. (2014). Exome resequencing reveals signatures of demographic and adaptive processes across the genome and range of black cottonwood (Populus trichocarpa). Mol. Ecol. 23, 2486–2499. doi: 10.1111/mec.12752
Zsuffa, L., Giordano, E., Pryor, L. D., and Stettler, R. F. (1996). “Trends in poplar culture: some global and regional perspectives,” in Biology of Populus and its Implications for Management and Conservation, eds R. F. Stettler, H. D. Bradshaw Jr., P. E. Heilman, and T. M. Hinckley (Ottawa, ON: NRC Research Press), 515–539.
Keywords: Populus trichocarpa, black cottonwood, common garden, local adaptation, regression tree
Citation: Oubida RW, Gantulga D, Zhang M, Zhou L, Bawa R and Holliday JA (2015) Partitioning of multivariate phenotypes using regression trees reveals complex patterns of adaptation to climate across the range of black cottonwood (Populus trichocarpa). Front. Plant Sci. 6:181. doi: 10.3389/fpls.2015.00181
Received: 28 May 2014; Accepted: 06 March 2015;
Published: 27 March 2015.
Edited by:Heikki Hänninen, University of Helsinki, Finland
Reviewed by:Glenn Thomas Howe, Oregon State University, USA
Outi Savolainen, University of Oulu, Finland
Copyright © 2015 Oubida, Gantulga, Zhang, Zhou, Bawa and Holliday. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jason A. Holliday, Department of Forest Resources and Environmental Conservation, Virginia Polytechnic Institute and State University, 304 Cheatham Hall, Blacksburg, VA 24061, USA email@example.com