Pinus strobiformis is a widely distributed white pine, ranging from central and southern Arizona, New Mexico, western Texas, and southwestern Colorado in the United States, to the Sierra Madre Occidental in Mexico (Figure 1; Looney and Waring, 2013; Villagómez-Loza et al., 2014). According to the Mexican National Forest Inventory (2004–2009), performed by the Mexican National Forestry Commission (CONAFOR, 2009), P. strobiformis occurs within an area of about 2.3 million hectares in Mexico, mainly at elevations of between 2,500 and 3,000 m, although the suitable habitat is estimated to constitute a much larger area (Aguirre-Gutiérrez et al., 2015; but see Shirk et al., 2018). P. strobiformis occurs at elevations between 1,800 and 3,400 m in the southwestern United States (Looney and Waring, 2013) but occupies less than 2 % of the total forested area in Arizona (O’Brien, 2002).

In 2015, we collected cones and seeds from 65 P. strobiformis sites, located across a variety of abiotic conditions and a wide latitudinal gradient in Mexico and the United States and representing the current geographic distribution (Figure 1, Supplementary Tables S1–S5). Sites were selected based on accessibility and availability of ripe cones and seeds, with a minimum distance of 5.1 km between sites. All stands were situated in closed forests with minimal human disturbance (e.g., roads, cattle grazing, agriculture). Three to five trees with ripe cones were sampled at each site (a total of 297 trees). In each site, the sampled trees were separated by a minimum distance of 50 m, to minimize the chance of sampling the same genetic family.

For each tree we recorded: Latitude, longitude and elevation (m), diameter at breast height (cm), total tree height (m), first live branch height (m), and crown length (m, obtained from the difference between total height and first live branch height) (Supplementary Tables S1, S2).

A minimum of 15 ripe cones with no signs of the presence of insects or diseases were collected from each tree, stored and transported in labeled bags (total 4,455 cones). These cones were allowed to dry at ambient temperature and humidity until they opened, and the seeds were then extracted. Empty seeds were separated from filled seeds, either manually or with a blower. The filled seeds were then weighed in 10-seed lots (g) on a balance (González-Ávalos et al., 2006; Iglesias et al., 2012).

Ten ripe cones were selected at random from the 15 cones (or more) harvested from each tree, for measurements. Total length (cm) and width (cm) of the widest part of each cone were recorded (the width was measured with the angle of the cone facing upward and taking care to avoid exerting pressure on the cone) (Bramlett et al., 1977; González-Ávalos et al., 2006; Iglesias et al., 2012; Prieto-Ruíz et al., 2014). The cone angle was measured with a 360° circular protractor, by placing the upper stem of the cone aligned with 0° (Forde, 1964; Wheeler and Guries, 1982).

Three of the ten cones per tree were randomly chosen for further measurements; three scales of the upper and three of the lower half of each cone were selected and removed. These scales were chosen from each third of the cone circumference, at angles of approximately 0°, 120°, and 240°. A total of 2,673 scales from the upper and 2,673 from the lower part of the cones were gotten from the 297 trees. The length and width of each scale was then measured with a digital caliper (mm). The length of the scale to the tip and the width at the widest part of the scale were recorded. The angle of each scale was measured using a circular protractor and aligning each scale to 90°. Cone specific trait means were calculated for each tree and site.

For each site (Figure 1), the following abiotic and biotic factors were recorded to examine their impact on the morphological traits of P. strobiformis: geographical aspect (0° to 360°), slope (%), occurrence of P. strobiformis regeneration (presence), presence or not of Ribes spp. (an alternative host for the white pine blister rust pathogen) and presence or not of woody species in the neighborhood (Supplementary Tables S4, S5). We only recorded the presence of these species within a 20 m radius of each individual P. strobiformis tree and calculated the relative frequency of these species. The sampled trees were also examined for signs and symptoms of white pine blister rust, and the following rating system was applied when the disease was present: cankers present on branches (1), bole (2), or bole and branches (3) (Arvanitis et al., 1984); these ratings were then transformed into frequency of occurrence prior to analyses (Supplementary Table S2). The climate data (20 temperature and precipitation variables) were downloaded from the PRISM database for the period from 1961 to 1990¹ (Supplementary Table S3). For further analysis, the geographical aspect (0° to 360°) was transformed into a cosine index – as 0° and 360° have the same zenithal aspect (Cushman and Wallin, 2002). Occurrence data were transformed into frequency of occurrence prior to analyses of the presence of trees and shrubs (Supplementary Tables S4, S5). We assumed that the presence of other tree species around P. strobiformis does not directly affect its morphological traits. However, this factor was used as additional proxy of other abiotic factors (such as soil traits) possibly influencing the morphological traits of P. strobiformis (Zhang et al., 2016).

Since P. strobiformis x P. flexilis hybrids, which were only reported in the US populations, could influence morphological traits, we also recorded the hybrid degree (the relative hybrid frequency or proportion of P. strobiformis to P. flexilis per stand) as predictor variable, following Menon et al. (2018) (Supplementary Table S1).

We implemented redundancy analysis with a multi-tiered variance partitioning method (e.g., Cushman and McGarigal, 2002), in the “vegan” R package (version 3.3.4) (Oksanen et al., 2010) to quantify the independent and joint ability of each set of environmental variables to predict the morphological characteristics of all sampled white pine individuals, and to measure the importance and joint interactive effects of the variables together. Importantly, this also enabled us to quantify the amount of morphological variation not explained by environmental variation in our data set. We computed a four-way partitioning among all climatic, vegetation, geographical and topographical variables. The spatial variables included the eigth spatial trend surface analysis terms (e.g., x, y, x^2, y^2, xy, x^2y, xy^2, and x^2y^2; Cushman and McGarigal, 2002).

Selection of appropriate independent variables is fundamental for achieving effective predictive models (Krzanowski, 1987; Gnanadesikan et al., 1995; Maronna et al., 2018). We tested for significant differences in the mean values of explanatory environmental variables across the 65 P. strobiformis stands with morphological data (Supplementary Tables S1–S5) to describe variation across populations. Many of our predictor variables were not normally distributed (the variables are listed in Supplementary Tables S1–S5). Due to the absence of normality, we used the non-parametric Kruskal-Wallis test (Kruskal and Wallis, 1952) to evaluate whether the observed median’s differences in independent variables between P. strobiformis populations were statistically significant. All environmental variables for which we detected significant differences in median values of morphological traits (α = 0.01) were included in further analysis.

Additionally, the most important environmental factors influencing morphological variation were also determined using partial least squares and area under the univariate receiver operator curve, using the “varImp” function applied to the results of univariate analysis with the machine learning algorithm, using Random Forest [“caret” package, function train, methods “rf,” implemented in R (version 3.3.4; R Development Core Team, 2017)]. In each case, the unsupervised correlation filter was applied to the predictors prior to modeling (see details in the section “Regression techniques”).

The variable-importance measure was determined using partial least squares (Kuhn, 2012), on the basis of the weighted sums of the absolute regression coefficients. The weights are a function of the reduction of the sums of squares across the number of partial least squares components and are computed separately for each outcome. The contributions of the coefficients are thus weighted proportionally to the reduction in the sums of squares. The trapezoidal rule was used to compute the area under the receiver operator curve, which was used as a measure of the variables importance (Kuhn, 2012).

Finally, we used the non-parametric Spearman’s coefficient (r_s) to determine correlations between the variables and to estimate collinearity between important independent variables (selected by the Kruskal-Wallis test, partial least squares or receiver operator curve). When the r_s absolute value for the difference between two variables was greater than 0.70, only the variable with the lowest p value in the Kruskal–Wallis test was included in the regression models (as reported by Salas et al., 2017 and Shirk et al., 2018). The relationships between the most important variables and two most spatially dependent cone and seed traits were represented graphically.

The number of variables was determined by the rule of ten events per variable (McGarigal et al., 2000; Vittinghoff and McCulloch, 2007, i.e., a maximum of the six most important variables for the 65 P. strobiformis stands were included in the models). Regression models, including 5-fold cross-validations, were used to predict the most spatially dependent P. strobiformis cone and seed traits for selected important, independent variables in each stand. Six machine learning algorithms were implemented in the “caret” package and function “train” models: (i) linear regression (method = “lm”), (ii) Random Forest (method = “rf”), (iii) Neural Network (method = “nnet”), (iv) Model Averaged Neural Network (method = “avNNet”), (v) Multi-Layer Perceptron (method = “mlpWeightDecay”), and (vi) Bayesian Regularized Neural Networks (method = “brnn”) (Venables and Ripley, 1999; Williams et al., 2018, http://topepo.github.io/caret/index.html) in R (version 3.3.4) (R Development Core Team, 2017). The goodness-of-fit of the regression model was evaluated by using the (pseudo) coefficient of determination (R²), root of the mean square error (RMSE), and mean squared error.

The values of all studied morphological traits increased from the northern (United States) to the southern (Mexican) populations (e.g., cone length and seed weight were larger in southern populations; Figure 2; Supplementary Figure S1). The best kriging model for cone length used Stein’s parameterization (R_k² = 0.89; URMSE = 1.72 cm; Supplementary Table S6). The goodness-of-fit values for the seed weight and scale top angle were slightly lower (R_k² = 0.75, URMSE = 0.04 g; R_k² = 0.74, URMSE = 3.04 mm). In general, the mean cone size, seed weight and scale angle were larger in the southern populations. The worst-performing spatial model was that for scale top width (R_k² = 0.05). There was a marginally positive relationship between latitude and DBH (R_k² = 0.06, p = 0.033), but there was no association between latitude and tree height (R_k² = 0.15, p = 0.38) (Supplementary Table S6).

Collectively, vegetation, spatial, topographic, geographical, and climatic variables explained 54.7% of the total variation in the morphological characteristics among sampled trees (Figure 3). Of the marginal effects of these groups of variables, geographic variables showed the best ability to predict morphological characteristics, accounting for 51.5% of the variance (Supplementary Table S7), followed by vegetation variables (43.8%), climate variables (43.4%) and finally topographic variation (36.9%). There was very high covariation in the explanation among these different sets of variables, with no set explaining more than 6% (set geographic) of the variation in morphological traits independent of the other variable sets. The climate and topographic variable sets had very close to zero (0.004) independent explanatory power. The largest component of explained variance (33.6%) was the four way interaction of all four variable sets, suggesting that there is strong co-variation in environmental, climate and geographical variables in their relationship to morphological characteristics of P. strobiformis across its range. The second largest variance component was the three-way interaction between climate, vegetation and geographic variable sets. This shows that in total 39.6% of the variance in morphological traits across the range of southwestern white pine are jointly predicted by simultaneous variation in climate, vegetation and geographic variable sets (Figure 3).

For the 65 stands with morphological data, our modeling results indicated that growing season precipitation (GSP) was the most important independent variable for predicting both cone length and seed weight (Supplementary Figures S2, S3). However, the independent variables that together produced the best model of cone length of P. strobiformis were the GSP, winter precipitation (WINP), summer precipitation balance (SMRPB), frequency of occurrence of several overstory tree species (Pseudotsuga menziesii (Mirbel) Franco, Pinus cooperi C.E. Blanco and Arbutus xalapensis Kunth), yielding an RMSE of 1.75 cm using the Random Forest (rf) (Table 1; Supplementary Figure S2). The second-best model for cone length (RMSE = 1.83 cm) was produced by the Bayesian Regularized Neural Networks (brnn) and the same independent variables.

In contrast, the GSP, frequency of occurrence of Pseudotsuga menziesii, mean annual temperature (MAT), SMRPB, frequency of occurrence of P. arizonica Engelm. and J. deppeana Steud. in the same site, were the variables that together provided the best prediction of seed weight of P. strobiformis [with an RMSE of 0.039 g using linear regression, (lm)]. The second best model of seed weight (RMSE = 0.040 g) was produced using the variables GSP, MAT, frequency of occurrence of Pseudotsuga menziesii, frequency of occurrence of P. strobiformis, P. arizonica and SMRPB, the brnn approach (Table 2; Supplementary Figures S3). Higher GSP and SMRPB, lower frequency of occurrence of Pseudotsuga menziesii, corresponded to longer mean cone length (Supplementary Table S8). Higher GSP, MAT, higher frequency of occurrence of P. arizonica and J. deppeana, and lower frequency of occurrence of P. menziesii and earlier SMRPB were correlated with greater mean seed weights (Supplementary Table S9).

The tree dimensions of DBH, total height and crown length were only marginally and non-significantly correlated with mean cone length (R² = 0.03, p = 0.15; R² = 0.02, p = 0.25; R² = 0.02, p = 0.25) and seed weight (R² = 0.07, p = 0.04; R² = 0.003, p = 0.67; R² = 0.03, p = 0.19), respectively. But, there was a moderately negative relationship between hybrid degree and both mean cone length and seed weight (R² = 0.30, p<0.00001), respectively.

Figure 4 shows the associations between mean cone length and seed weight and the most important variables: growing season precipitation (GSP), frequency of occurrence of P. menziesii in the neighborhood, and summer precipitation balance (SMRPB).

The regression results also illustrated how environmental (GSP, SMRPB, WINP, and MAT) factors influenced the cone length and seed weight of P. strobiformis. The largest cones and heaviest seeds were found in more humid and temperate climates, which are located in the southern half of the Mexican Sierra Madre Occidental. This result is in contrast to earlier work in other tree species in the United States by Baker (1972); Schimpf (1977), Sorensen and Miles (1978) and Stromberg and Patten (1990) who found that seeds were larger in drier sites. This difference may reflect the much larger and broader sample of our study, which covered the entire species range and therefore provided a more complete picture of phenotypic variation across the entire distribution than was possible with more limited prior studies. Conversely, Mazer (1989) and Telenius and Torstensson (1991) found no significant association between seed weight and moisture availability; again, those studies were limited by the extent and size of sampling. Leishman and Westoby (1994b) reported that the comparative evidence of an association between large seeds and dry habitats is very limited, despite the general assumption made in the scientific literature. Our study resolves this issue by evaluating a larger, more representative and range-wide sample, and shows strong associations of seed weight and cone length with latitude and climate, with larger seeds and longer cones in the southern, wetter part of the species’ range. However, hybridization with P. flexilis also leads to smaller cones and seeds of P. strobiformis at its northern border (Menon et al., 2018). Finally, hybridization effects with P. ayacahuite should not result in bigger seeds of P. strobiformis at the southern distribution border because seeds from P. ayacahuite are much smaller. However, hybridization with P. ayacahuite may lead to longer cones (Leal-Sáenz et al., 2020).

The strong covariation between topographical, geographic, climatic, and vegetation in their relationship to morphological traits suggests that the major variation in the measured morphological traits is associated with the simultaneous and covarying influences of spatial, climatic and topographic factors. The dominant explanatory power of climate, vegetation and geographic variable sets suggest that these factors are particularly important in their relationship to morphological characteristics, while variables in the topographic group have relatively weaker association with morphological traits. The inability of this large empirical sample to statistically separate the effects of climate, vegetation and geographic groups suggests that it is difficult to identify which specific environmental, spatial or climatic factors may be driving observed morphological variation.

The results reported here, however, are useful in showing that a majority of the variability in P. strobiformis morphological traits is explainable by climatic, topographical, spatial and vegetation variables and the degree of hybridization with P. flexilis. There are strong patterns of morphological variation in this species that are largely explainable by the joint and simultaneous effects of multiple factors. It is not surprising that there is high and inseparable covariation among these factors, as climatic variables are known to change in predictable ways with geographical location (latitude) and topography (elevation). Thus, one hypothesis to explore in common garden experiments is that the variation in morphological traits across the species range are primarily local adaptations to climate, which covaries with elevation and other ecological variables, and that the high apparent relationship between morphology and topography, space and vegetation community are spurious correlations with the actual climatic drivers. A second, alternative, hypothesis could be that competition with other tree species drives the morphological differences across the range, and that the distribution of these competitors is associated with climatic gradients, leading to covariation of morphological traits with climate. Our evidence suggests that cone and seed morphology align with increasing hybridization with P. flexilis, resulting in shorter cones and smaller seeds (Steinhoff and Andresen, 1971). Hybridization occurs primarily in the United States populations to the north, creating complex interactions between morphology, climate and genetics (Menon et al., 2018).

Various studies have reported that seeds are often larger under competitive conditions due to the enhanced survival of larger seeds (Salisbury, 1942, 1974; Grime and Jeffrey, 1965; Hutchinson, 1967; Ng, 1978; Foster and Janson, 1985; Mazer, 1989; Leishman and Westoby, 1994a). Egli (2017) argued that variation in seed weight is mainly related to variation in the rate of seed growth during rapid seed filling. Decades of provenance trials in forest trees provide evidence for wide variation in several key ecological traits (Alberto et al., 2013; Lind et al., 2018). For instance, in whitebark pine (Pinus albicaulis), Bower and Aitken (2008) suggested that the phenotypic variation is due to genetic and geographic differentiation that reflects the long-term adaptive evolution from the last glacial maximum and from local environmental adaptation. P. strobiformis has a large north to south range that results in a gradient of increasing summer precipitation and temperature of the coldest month. The latitudinal climatic gradient along with other abiotic factors also influenced the presence of other tree species around P. strobiformis. Therefore, the frequency of occurrence of these species corresponded to the cone length and seed weight as proxy of other abiotic and ecological factors of morphological traits (Zhang et al., 2016, Bañares-de-Dios et al., 2020). For example, we detected Pseudotsuga menziesii only in the United States sites, while Juniperus deppeana, P. arizonica, P. cooperi, and A. xalapensis were found only in the Mexico sites, and in the seed weight model the former species is associated with lower seed weights and the latter species with heavier average seed weights.

Many researchers have argued that variation in cone or seed size within species are closely associated with fitness and, therefore, with adaptive evolution (Roach, 1987; Winn, 1988; Biere, 1991; Platenkamp and Shaw, 1993; Ji et al., 2011). Seed size may also influence seed dispersal (Herrera et al., 1994; Jordano, 1995; Martínez et al., 2007; Shimada et al., 2015). Wang and Ives (2017) reported that seed weight affected almost all choices that rodents made in eating, removing and storing individual seeds. At the level of individual trees, larger seeds have improved probabilities of both predation and effective dispersal. Leslie et al. (2017) also reported that larger seeds were dispersed by animals while smaller seeds were dispersed by the wind and that the morphology of the cones is associated with the size of the seeds. Other studies have reported that latitude, genome size, forest structure, growth form and seed dispersal are related to differences in seed size (Salisbury, 1974; Lord et al., 1997; Leishman et al., 2000; Beaulieu et al., 2007). As the climate continues to shift, the morphological differences between the northern and southern P. strobiformis populations may result in differing reforestation patterns and affect management recommendations (Goodrich et al., 2018; Schoettle et al., 2018).

Experimental common garden studies could separate these influences to some degree, by replicating species combinations across climatic gradients. Common garden experiments are important additionally to quantify the relative degree of phenotypic plasticity and, while they are unlikely to show local adaptation in action, they can show the signature of past local adaptation in the portion of the variance that is not ascribable to phenotypic plasticity. Ultimately, the observed morphological variation is a product of demographic history, local adaptation to climate and environment, phenotypic plasticity, and how these interact in the context of local environmental conditions. The importance of demographic history and genetic drift could be noted in the high explanatory power of geographic variables even though several of them exhibited co-variation with climate and vegetation. To fully resolve the drivers of morphological variation, therefore, replicated and controlled common-garden experiments across the species range, representing combinations of all variable sets (e.g., topography, space, environment and genetics) are essential (Cushman, 2014). These common garden experiments should reciprocally transplant tree families from across the full range, and associate their expressed phenotypic variation with garden specific values of environmental variables, the values of these variables in the location of their maternal trees, and their individual genomic characteristics. Only through replicated and controlled experiments across ecological gradients that associate expressed variation with both environmental and genomic variability can we reliably identify and quantify the drivers of phenotypic variation (Sork et al., 2013).

In practice, however, such common garden studies may have limited practicality given the length of time needed for common garden plots to establish and mature to such a degree that they influence the expression or selection of morphological traits of long lived tree species. In addition, these gardens would have to be large to sufficiently represent forest stand conditions that would reflect competitive dynamics that drive selection or phenotypic plasticity. Finally, given that potentially a large amount of the variation in expressed phenotype is genetically controlled local adaptation, such common gardens would not necessarily show changes in phenotype in time frames feasible for experimentation, given multiple generations of selection are required to observe microevolutionary change in response to selection pressures related to competition or climate. Maintaining a large sample of spatially extensive gardens for many decades is a logistical and financial challenge.

Broad-scale sampling of genomic variation, morphological characteristics, and environmental variables across the species range provides another powerful framework to disentangle space, environment, climate, competition and genetic factors in influencing morphology (e.g., Cushman, 2014). Therefore, it is likely that the most effective way to separate the covarying influences of climate and competition on phenotypic variation will be through large-scale sampling of genomic variation across the population and association of this genomic variation with local environmental, climatic, and community conditions, and partitioning of the variance in morphological traits that are explainable by genomic variation vs. environmental variation. This would show the portion of morphological variation that is genetically correlated, which potentially indicates the degree of local adaptation. Specifically, while sampling standing variation in morphology and genomics across a species range does not enable experimental control necessary to isolate particular drivers with strong inferences, it does allow comparative mensurative designs that can trade space for time, and, critically, sample conditions in situ and at scales in both space and time that are operative and influential on the evolution and community dynamics of trees (McGarigal and Cushman, 2002). Thus we recommend future research combine and couple both broad-scale large-sample analysis of gene-environment gradients across the species range that can describe variation at broad scales and over long timer periods, with targeted, replicated and controlled common garden experiments that can isolate particular drivers of variation at small scales and over short times (Cushman, 2014).

Given the difficulty of reliably isolating drivers and apportioning explanatory variance among them with either observational or common garden experiments on long-lived tree species, we also recommend the use of simulation modeling (e.g., Landguth and Cushman, 2010; Landguth et al., 2017). Specifically, employing an individual-based, spatially explicit eco-evolutionary model (e.g., Landguth et al., 2020) to simulate the interactions of different degrees of gene flow, drift, environmental selection and phenotypic plasticity provides a unique means to explore the potential interactions of these factors and quantify the patterns of genomic and phenotypic variation that can be expected under these interactions (e.g., Cushman, 2015; Cushman and Landguth, 2016).

The ability of P. strobiformis to colonize its expected future range (Shirk et al., 2018) will be influenced by numerous factors, including colonization at the leading edge of the range shift, seed dispersal dynamics, resistance to white pine blister rust, competition with other species, introgression and hybridization (Menon et al., 2020), and genetic adaptation to local climate (Goodrich et al., 2016; Bucholz et al., 2020). Through our results, we have a better understanding of the environmental controls on cone and seed morphology, and we can more adequately evaluate seed provenances and transfer zones and provide better information for assisted migration strategies. Current best practices for seed transfer, such as (a) promoting a tight network of seed stands to prevent greater loss of local genetic variants and structure and (b) using the seeds to establish seedlings within a limited radius from each seed stand/provenance (Hernández-Velasco et al., 2017), are likely to exhibit continued success in the near future. Such practices are essential for current and near-future reforestation programs, including assisted migration or the establishment of new populations in areas that should be appropriate for specific species under expected climate change scenarios (Wehenkel et al., 2017). However, more innovative practices will be required under new climates not conducive to the reproductive success of local populations. This is particularly important given the large projected climate-driven range shift of P. strobiformis and predicted large climatic changes within the parts of the range that are likely to remain occupied (Shirk et al., 2018). We recommend a proactive approach, in which reforestation programs incorporate both local and regional seed sources and allow for assisted gene flow (Aitken and Bemmels, 2016). For example, a manager of P. strobiformis might include seed sources from the provenance that is most similar to the projected future climate of the planting location (given climate predictions for the species; Shirk et al., 2018) alongside local sources on sites of least aridity (Bucholz et al., 2020). Such a strategy promotes the local genetic structure while providing for adaptation to changing climates and allowing managers flexibility in determining seed sources most suitable given other management considerations (e.g., white pine blister rust). As cone length and seed weight are probably also related to local/regional genetic adaptation (Rawat and Bakshi, 2011), these variables could be included in future climate-based seed transfer guidelines, providing suggestions for average cone lengths and seed weights most appropriate for a given site, including information on degree of hybridization present or desired for a planting site. Given the observed high covariation between climatic, community structure and geographic factors in this study, and the difficulty in separating them statistically or experimentally, there may be no silver bullet for managers to decide on the merits of different planting and assisted migration strategies. Again, given the difficulty of using observational or experimental data to guide these decisions, it may be useful to augment such studies with simulations that can control the interactions of these factors across large ranges of scale in both space and time to quantify the potential roles and influences of each factor in the context and interactions with the others (e.g., Landguth et al., 2017).

Importantly, these results show that the pattern of morphological variation in P. strobiformis is highly predictable in relation to a combination of geographic, floristic, climatic and topographical variables. It is interesting to note that our analysis did not include several environmental variables that are known to strongly predict the distribution of P. strobiformis (e.g., soils; Shirk et al., 2018). Thus, it is likely that including a broader set of environmental variables that are limiting to P. strobiformis fitness at different scales and in different contexts would increase the amount of variance explained in morphological characteristics across the species range. Additionally, our analysis does not formally integrate observed genomic variation among individual trees in comparison to expressed phenotype across gradients of environmental, climatic and geographic gradients. Future work should therefore focus on collecting these factors simultaneously for a large number of trees across the full extent of the species range and ecological conditions to enable more rigorous evaluation of the degree to which genomic variation can explain observed phenotypic variation, and to what degree it is covarying and thus potentially controlled by selection along environmental gradients. Such studies would also be able to quantify the degree to which admixture with peripatric sister species (e.g., P. flexilis) confounds and contributes to observed variation in phenotypic characteristics along geographical and environmental gradients. Another area that would be valuable to integrate into analyses of relationships between phenotypic variation, geography, community structure and environmental gradients would be formal accounting of the influences of phylogeographic and demographic history. Phylogeographic and demographic history create non-stationary and non-equilibrium patterns of genetic structure across populations that are not linearly related to local adaptation or patterns of gene flow (e.g., Dyer et al., 2010). Developing and integrating methods that can account for this therefore is valuable. Again, simulation modeling (Cushman, 2015) may be the most powerful way to account for the interactions of gene flow, drift and selection within varying contexts of phylogeographic and demographic history, given their ability to stipulate and control all the factors that interact in a way that enables simulation experiments to robustly evaluate each factor and its interactions with others.