Morphological Differences in Pinus strobiformis Across Latitudinal and Elevational Gradients

Abstract

The phenotype of trees is determined by the relationships and interactions among genetic and environmental influences. Understanding the patterns and processes that are responsible for phenotypic variation is facilitated by studying the relationships between phenotype and the environment among many individuals across broad ecological and climatic gradients. We used Pinus strobiformis, which has a wide latitudinal distribution, as a model species to: (a) estimate the relative importance of different environmental factors in predicting these morphological traits and (b) characterize the spatial patterns of standing phenotypic variation of cone and seed traits across the species’ range. A large portion of the total variation in morphological characteristics was explained by ecological, climatic and geographical variables (54.7% collectively). The three climate, vegetation and geographical variable groups, each had similar total ability to explain morphological variation (43.4%, 43.8%, 51.5%, respectively), while the topographical variable group had somewhat lower total explanatory power (36.9%). The largest component of explained variance (33.6%) was the four-way interaction of all variable sets, suggesting that there is strong covariation in environmental, climate and geographical variables in their relationship to morphological traits of southwest white pine across its range. The regression results showed that populations in more humid and warmer climates expressed greater cone length and seed size. This may in part facilitate populations of P. strobiformis in warmer and wetter portions of its range growing in dense, shady forest stands, because larger seeds provide greater resources to germinants at the time of germination. Our models provide accurate predictions of morphological traits and important insights regarding the factors that contribute to their expression. Our results indicate that managers should be conservative during reforestation efforts to ensure match between ecotypic variation in seed source populations. However, we also note that given projected large range shift due to climate change, managers will have to balance the match between current ecotypic variation and expected range shift and changes in local adaptive optima under future climate conditions.

Introduction

The interest of evolutionary ecologists has long been focused on patterns of phenotypic variation along environmental gradients (e.g., latitude, altitude and climate). However, this variation is also determined by the covariance between genetic and environmental influences (Conover et al., 2009). Studying patterns of phenotypic variation among many individuals along broad ecological and climatic gradients is a powerful framework to understand the patterns and processes that govern phenotypic expression and variation (Endler, 1986). Early ecologists frequently noted that phenotypes of many species change predictably along large-scale gradients of latitude, altitude and water depth, providing the basis for several so-called ecological rules (Bergmann, 1847; Atkinson and Sibly, 1997). Morphological traits are important characteristics (Violle et al., 2007) that help distinguish entities at all levels of biological organization, including life cycles, ecological and geographical distributions, and evolutionary and conservation status (Kaplan, 2001; Gregorius et al., 2003). For example, Delagrange et al. (2004) showed that the crown morphology of Acer saccharum and Betula alleghaniensis is influenced by light availability and tree size. Wahid et al. (2006) found that the morphology of Pinus pinaster cones and needles vary predictably along both latitudinal and altitudinal gradients and Gil et al. (2002) reported elevation-dependent cone and seed variation in Pinus canariensis.

Forest trees are often adapted to environmental gradients at multiple spatial scales (Morgenstern, 1996; Savolainen et al., 2007). The phenotype of an individual is influenced by their genotype, their environment and interactions between them (de Jong, 1990; Falconer and Mackay, 1996). The evolutionary response of a phenotypic trait to selection depends on genetic control of the trait, heritability of the trait and differential fitness of different morphotypes of the trait in different environmental conditions (Price, 1970). Quantifying morphological variation within a widely distributed species across a large spatial extent is often necessary to identify the relative importance of different environmental factors in relation to variation in intraspecific morphological traits (Boyd, 2002; Herrera et al., 2002; Franks et al., 2014; Ji et al., 2016). Such knowledge can then be applied to climate change adaptation, through strategies such as assisted migration and assisted gene flow (e.g., Aitken and Bemmels, 2016).

White pine species (subgenus Strobus, section Quinquefoliae, subsection Strobus (Gernandt et al., 2005)) are widespread in the temperate forest ecosystems of North America (Kral, 1993; Gernandt and Pérez-de la Rosa, 2014), making them a suitable focal group to investigate how the phenotype is influenced by environmental gradients. Members of subgenus Strobus, section Quinquefoliae, subsection Strobus have five needles per bundle, such that species in this group are commonly referred to as five-needle pines. Several white pines are also characterized by large cones (Farjon and Styles, 1997) and large seeds (Frankis, 2009). These include the southwestern white pine (Pinus strobiformis Engelm.), which has wide variation in cone length, 20–50 cm (Frankis, 2009), seed length, 1.2–1.8 cm (Farjon and Styles, 1997), seed weight of 0.140–0.411 g (Leal-Sáenz et al., 2020), and a geographical range that includes the Sierra Madre Occidental in Mexico and the southwestern United States (Figure 1). Southwestern white pine is of both commercial and ecological value. Commercial products include timber, cellulose, resin and pulp (Farjon and Styles, 1997; Villalobos-Arámbula et al., 2014).

FIGURE 1

Ecologically, P. strobiformis is an important tree species in the ecosystems it inhabits (Moreno-Letelier and Piñero, 2009; Bower et al., 2011; Looney and Waring, 2013). The large and nutritious seeds of white pines frequently form the foundation of large trophic networks and promote enhanced biodiversity (Keane and Cushman, 2018). Seeds are an important source of food for many birds and mammals, including corvids, parrots, mice, voles, chipmunks, squirrels and bears (Samano and Tomback, 2003; Tomback and Achuff, 2010; Looney and Waring, 2013). In addition to contributing to biodiversity, P. strobiformis also provides watershed protection and serves as a major floristic component of montane forest ecosystems (Wehenkel et al., 2014), and is a key element of mixed conifer forests across its range (Looney and Waring, 2013; Shirk et al., 2018).

Pinus strobiformis is highly susceptible to white pine blister rust (WPBR), an invasive tree disease caused by the fungus Cronartium ribicola (Geils and Vogler, 2011), as are all other white pine species (Hoff and Hagle, 1990; Sniezko et al., 2008). The dual threats of this non-native fungal pathogen and the projected warmer, drier climate have created an uncertain future for P. strobiformis, as well as other white pine species (Landguth et al., 2017; Keane and Cushman, 2018; Shirk et al., 2018). A contraction of more than 60% in the species range by 2080 is predicted under some scenarios of greenhouse gas emissions (Shirk et al., 2018), whereas other scenarios include a northerly shift of more than 1,000 km in the mean latitude and an increase of 500 m in the mean elevation of suitable habitat for the species.

Robust measurement of the degree to which morphological variation is determined by non-genetically controlled developmental responses to environmental gradients requires a statistical analysis that can jointly and simultaneously account for the ability of environmental factors to predict morphological characteristics. RDA is a multivariate direct ordination approach that is ideally suited to predicting multiple response variables as a function of multiple predictor variables and is the most well-known and trusted statistical method to conduct multivariate variance partitioning (Van Den Wollenberg, 1977; McGarigal et al., 2000). By implementing a series of partial redundancy analysis ordinations, the independent and joint effects of multiple sets of predictor variables can be separated (e.g., Borcard et al., 1992; Cushman and McGarigal, 2002, 2004). For instance, the proportion of morphological trait variance explained by pure effects of climatic data, spatial location or topography and the interaction between variables can be quantified. A series of partial redundancy analyses can then be used to identify environmental gradients correlated with phenotypic variation across the species range. This partitioning of the variance in morphological traits that can be explained independently and jointly by climatic and geographic factors is essential to understand the factors that drive morphological variation and the scales at which they operate.

Faced with an increasingly complex and rapidly changing future, forest managers need science-driven strategies to maintain tree species and forests. Given the ecological importance of Pinus strobiformis and the multiple challenges to future regeneration of the species, it is critical to understand factors underlying the morphological trait variation of reproductive structures (e.g., cones and seeds) in order to better manage reforestation and restoration efforts in the future. We used extensive sampling across the full range of Pinus strobiformis, to (a) estimate the relative importance of the effects of different environmental factors on these morphological traits and (b) characterize the spatial patterns of standing phenotypic variation of cone and seed morphological traits.

Materials and Methods

Species and Study Sites

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.

Tree, Cone and Seed Morphology Data Collection

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.

Environmental Factors and Hybrid Degree

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).

Statistical Analyses

Detecting Spatial Dependence in Morphological Traits by Ordinary Kriging Analysis

We conducted ordinary kriging (ordinary Gaussian process regression model) to spatially interpolate morphological traits, including a 10-fold point-by-point cross validation (Batista et al., 2016). The statistical software R (version 3.3.4) (R Development Core Team, 2017) and the Interpolation Kriging package (ArcGIS Desktop 10.5, 2016) were used to describe first-order variation in the spatial pattern of the morphological cone and seed traits, as well as the dasometric traits of the P. strobiformis sampled trees under study. We tested the following mathematical models for the semivariance: the spherical model, exponential model, Gaussian model and the Stein’s parameterization. The corrected coefficient of determination between the observed and predicted values (R²), the Unbiased Root Mean Squared Error of the residuals (URMSE), the mean squared error and the mean absolute error were used to assess the goodness-of-fit. Finally, the cone and seed traits with the best kriging model were selected for further regression models of morphological traits with respect to environmental variables. The modeling was carried out using the “SP” (Pebesma and Bivand, 2005) and “automap” packages (Hiemstra et al., 2009) including the CRS, SpatialPixelsDataFrame, autoKrige, autoKrige.cv, and compare.cv functions in R (version 3.3.4) (R Development Core Team, 2017).

Redundancy Analysis and Variance Partitioning

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 Independent Environmental Factors Influencing Morphological Traits

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.

Modeling Spatially Dependent Cone and Seed Traits by Machine Learning Regression Methods

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.

Results

Detecting Spatial Dependence in Morphological Traits by Ordinary Kriging Analysis

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).

FIGURE 2

Redundancy Analysis and Variance Partitioning

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).

FIGURE 3

Modeling Spatially Dependent Cone and Seed Traits Using Machine Learning Methods

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).

Associations of Cone Length and Seed Weight With Tree Dimension and Hybrid Degree

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).

FIGURE 4

Discussion

Modeling Spatially Dependent Cone and Seed Traits by Machine Learning Methods

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).

Redundancy Analysis and Variance Partitioning

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).

Common Gardens, Gradient Modeling and Simulation

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).

Implications for Management and Conservation

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).

Scope and Limitations

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.

Conclusion

We showed that most variation in P. strobiformis morphological characteristics is strongly correlated with climatic gradients, suggesting selection for different morphological characteristics under different climatic contexts. However, we could not determine how much of the morphological variation is driven by climate independently of covarying topographical, ecological and spatial factors. In addition, our study does not quantify the degree to which observed morphological variation is genetically controlled, nor how much phenotypic plasticity there is. We advocate for replicated and controlled reciprocal transplant common garden experiments (e.g., Sork et al., 2013; Cushman, 2014) which can enable rigorous separation of the amount of phenotypic variance controlled by genotype, by the environment and by the interaction between genotype and environment. In addition, we suggest coupling such common garden studies to broad-scale gradient modeling of genomic and morphological variation across geographic, climatic and community gradients (Cushman, 2014). Simulation modeling, which can control all the potential factors driving covariation between genomic, phenotypic, geographical, and environmental factors at scales in space and time relevant to population responses are likely to be critical to rigorously quantify these relationships and untangle them. Ultimately a combination of common garden, gradient modeling and simulation studies provides the best means for advancing the important and challenge task of understanding and predicting eco-evolutionary dynamics across broad populations in complex and dynamic environments (Cushman, 2014). The results presented here, however, are useful and important in showing spatial range-wide patterns of phenotypic variation that are strongly associated with environmental gradients, and that a large portion of this variation is also associated with the joint effects of climate, topography, latitude and regional vegetation community composition.

The combination of variance partitioning and machine learning algorithms implemented here provides a clear demonstration of both the relative importance and independent effects of different climatic, geographic and environmental factors in driving morphological variation in P. strobiformis, but also identify the main patterns of this variation and the variables that are most strongly associated with them. We find that cone length, seed weight and cone morphology are strongly related to temperature and precipitation and increase in warmer and wetter parts of the species range. These results show that spatial modeling across a species range can yield accurate predictions of morphological traits as a function of environmental gradients. These models predict the considerable differences in geographical, topographical, climate, adaptive and morphological variables in the species range and may help to distinguish the actual seed provenance of P. strobiformis.

References

• 1

  Aguirre-GutiérrezJ.Serna-ChavezH. M.Villalobos-ArambulaA. R.Pérez de la RosaJ. A.RaesN. (2015). Similar but not equivalent: ecological niche comparison across closely related Mexican white pines.Divers. Distribut.21245–257. 10.1111/ddi.12268

  • CrossRef
  • Google Scholar

• 2

  AitkenS. N.BemmelsJ. B. (2016). Time to get moving: assisted gene flow of forest trees.Evolut. Appl.9271–290. 10.1111/eva.12293

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  AlbertoF. J.AitkenS. N.AlíaR.González-MartínezS. C.HänninenH.KremerA.et al (2013). Potential for evolutionary responses to climate change–evidence from tree populations.Glob. Change Biol.191645–1661. 10.1111/gcb.12181

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  ArcGIS Desktop 10.5 (2016). © 1999-2016. Version: 10.5.0.6491. Property of Environmental Systems Research Institute, Inc. (Esri) Inc., All Rights Reserved. Redlands, California, United States.

  • Google Scholar

• 5

  ArvanitisL. G.GodbeeJ. F.PortaI. (1984). Pitch canker impact on volume growth: a case study in slash pine plantations.Southern J. Appl. For.843–47. 10.1093/sjaf/8.1.43

  • CrossRef
  • Google Scholar

• 6

  AtkinsonD.SiblyR. M. (1997). Why are organisms usually bigger in colder environments? Making sense of a life history puzzle.Trends Ecol. Evol.12235–239. 10.1016/S0169-5347(97)01058-6

  • CrossRef
  • Google Scholar

• 7

  BakerH. G. (1972). Seed weight in relation to environmental conditions in California.Ecology53997–1010. 10.2307/1935413

  • CrossRef
  • Google Scholar

• 8

  Bañares-de-DiosG.MacíaM. J.Granzow-de la CerdaÍArnelasI.de CarvalhoG. M.EspinosaC. I.et al (2020). Linking patterns and processes of tree community assembly across spatial scales in tropical montane forests.Ecology1011–13. 10.1002/ecy.3058

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  BatistaA. P. B.MelloJ. M. D.RaimundoM. R.ScolforoH. F.ReisA. A. D.ScolforoJ. R. S. (2016). Species richness and diversity in shrub savanna using ordinary kriging.Pesquisa Agropecu. Brasil.51958–966. 10.1590/S0100-204X2016000800008

  • CrossRef
  • Google Scholar

• 10

  BeaulieuJ. M.MolesA. T.LeitchI. J.BennettM. D.DickieJ. B.KnightC. A. (2007). Correlated evolution of genome size and seed mass.New Phytol.173422–437. 10.1111/j.1469-8137.2006.01919.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  BergmannC. (1847). Über die Verhältnisse der Wärmeökonomie der Thiere zu ihrer Größe.Gött. Stud.1595–708.

  • Google Scholar

• 12

  BiereA. (1991). Parental effects in Lychnis floscuculi. I: seed size, germination and seedling performance in a controlled environment.J. Evol. Biol.4447–465. 10.1046/j.1420-9101.1991.4030447.x

  • CrossRef
  • Google Scholar

• 13

  BorcardD.LegendreP.DrapeauP. (1992). Partialling out the spatial component of ecological variation.Ecology731045–1055. 10.2307/1940179

  • CrossRef
  • Google Scholar

• 14

  BowerA. D.AitkenS. N. (2008). Ecological genetics and seed transfer guidelines for Pinus albicaulis (Pinaceae).Am. J. Bot.9566–76. 10.3732/ajb.95.1.66

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  BowerA. D.McLaneS. C.EckertA.JorgensenS.SchoettleA.AitkenS. (2011). “Conservation genetics of high elevation five-needle white pines,” in Proceedings of the High Five Symposium and RMRS-P-63 on The Future of High-Elevation, Five-Needle White Pines in Western North America, edsKeaneR. E.TombackD. F.MurrayM. P.SmithC. M. (Fort Collins, CO: U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station), 98–117.

  • Google Scholar

• 16

  BoydA. (2002). Morphological analysis of Sky island populations of Macromeria viridiflora (Boraginaceae).Syst. Bot.27116–126. 10.1043/0363-6445-27.1.116

  • CrossRef
  • Google Scholar

• 17

  BramlettD. L.BelcherE. W.Jr.DebarrG. L.HertelG. D.KarrfaltR. P.LantzC. W.et al (1977). Cone Analysis of Southern Pines- A Guidebook. General Technical Report SE 13.Asheville, NC: U.S. Department of Agriculture, Forest Service, Southeastern Forest Experiment Station.

  • Google Scholar

• 18

  BucholzE.WaringK. M.KolbT. E.SwensonJ.WhippleA. (2020). Water relations and drought response of Pinus strobiformis.Can. J. For. Res.50905–916. 10.1139/cjfr-2019-0423

  • CrossRef
  • Google Scholar

• 19

  CONAFOR (2009). Comisión Nacional Forestal.México: National Forestry Commission of Mexico.

  • Google Scholar

• 20

  ConoverD. O.DuffyT. A.HiceL. A. (2009). The covariance between genetic and environmental influences across ecological gradients.Ann. N. Y. Acad. Sci.1168100–129. 10.1111/j.1749-6632.2009.04575.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  CushmanS. A. (2014). Grand challenges in evolutionary and population genetics: the importance of integrating epigenetics, genomics, modeling, and experimentation.Front. Genet.5:197. 10.3389/fgene.2014.00197

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  CushmanS. A. (2015). Pushing the envelope in genetic analysis of species invasion.Mol. Ecol.24259–262. 10.1111/mec.13043

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  CushmanS. A.LandguthE. L. (2016). Spatially heterogeneous environmental selection strengthens evolution of reproductively isolated populations in a Dobzhansky–muller system of hybrid incompatibility.Front. Genet.7:209. 10.3389/fgene.2016.00209

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  CushmanS. A.McGarigalK. (2002). Hierarchical, multi-scale decomposition of species-environment relationships.Landscape Ecol.17637–646. 10.1023/A:1021571603605

  • CrossRef
  • Google Scholar

• 25

  CushmanS. A.McGarigalK. (2004). Patterns in the species–environment relationship depend on both scale and choice of response variables.Oikos105117–124. 10.1111/j.0030-1299.2004.12524.x

  • CrossRef
  • Google Scholar

• 26

  CushmanS. A.WallinD. O. (2002). Separating the effects of environmental, spatial and disturbance factors on forest community structure in the Russian Far East.For. Ecol. Manag.168201–215. 10.1016/S0378-1127(01)00744-7

  • CrossRef
  • Google Scholar

• 27

  de JongG. (1990). Genotype-by-environment interaction and the genetic covariance between environments: multilocus genetics.Genetica81171–177. 10.1007/BF00360862

  • CrossRef
  • Google Scholar

• 28

  DelagrangeS.MessierC.LechowiczM. J.DizengremelP. (2004). Physiological, morphological and allocational plasticity in understory deciduous trees: importance of plant size and light availability.Tree Physiol.24775–784. 10.1093/treephys/24.7.775

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  DyerR. J.NasonJ. D.GarrickR. C. (2010). Landscape modelling of gene flow: improved power using conditional genetic distance derived from the topology of population networks.Mol. Ecol.193746–3759. 10.1111/j.1365-294X.2010.04748.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  EgliD. B. (2017). Seed Biology and the Yield of Grain Crops, 2nd Edn. Lexington, KY: University of Kentuchy.

  • Google Scholar

• 31

  EndlerJ. A. (1986). Natural Selection in the Wild. Monographs in Population Biology 21.Princeton, NJ: Princeton University Press.

  • Google Scholar

• 32

  FalconerD. S.MackayT. F. C. (1996). Introduction to quantitative genetics, 4th Edn. London: Pearson Education Ltd.

  • Google Scholar

• 33

  FarjonA.StylesB. T. (1997). Pinus (pinaceae).Flora Neotropica. Monograph 75. Bronx, NY: The New York Botanical Garden.

  • Google Scholar

• 34

  FordeM. B. (1964). Variation in natural populations of Pinus radiata in California, New Zealand.J. Bot.2459–485. 10.1080/0028825X.1964.10428765

  • CrossRef
  • Google Scholar

• 35

  FosterS.JansonC. H. (1985). The relationship between seed size and establishment conditions in tropical woody plants.Ecology66773–780. 10.2307/1940538

  • CrossRef
  • Google Scholar

• 36

  FrankisM. (2009). The High Altitude White Pines of Mexico and the Adjacent SW USA. (Pinus L. Subgenus Strobus Lemmon, Pinaceae). Available online at: http://pinetum.org/articles/Pinus/PNstylesiiFrankis2009.pdf(accessed August 19, 2019).

  • Google Scholar

• 37

  FranksS. J.WeberJ. J.AitkenS. N. (2014). Evolutionary and plastic responses to climate change in terrestrial plant populations.Evolut. Appl.7123–139. 10.1111/eva.12112

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  GeilsB. W.VoglerD. R. (2011). “A natural history of Cronartium ribicola,” in The Future of High-Elevation, Five-Needle White Pines in Western North America, edsKeaneR. E.TombackD. F.MurrayM. P.SmithC. M. (Fort Collins, CO: US Department of Agriculture, Forest Service, Rocky Mountain Research Station), 210–217.

  • Google Scholar

• 39

  GernandtD. S.LópezG. G.GarcíaS. O.ListonA. (2005). Phylogeny and classification of Pinus.Taxonomy5429–42.

  • Google Scholar

• 40

  GernandtD. S.Pérez-de la RosaJ. A. (2014). Biodiversidad de Pinophyta (coníferas) en México.Rev. Mexic. Biodiver.85126–133. 10.7550/rmb.32195

  • CrossRef
  • Google Scholar

• 41

  GilL.ClimentJ.NanosN.MutkeS.OrtizI.SchillerG. (2002). Cone morphology variation in Pinus canariensis Sm.Plant Syst. Evol.23535–51. 10.1007/s00606-002-0218-9

  • CrossRef
  • Google Scholar

• 42

  GnanadesikanR.KettenringJ. R.TsaoS. L. (1995). Weighting and selection of variables for cluster analysis.J. Classific.12113–136. 10.1007/BF01202271

  • CrossRef
  • Google Scholar

• 43

  González-ÁvalosJ.García-MoyaE.Vargas-HernándezJ. J.Trinidad-SantosA.Romero-ManzanaresA.Cetina-AlcaláV. M. (2006). Evaluación de la producción y análisis de conos y semillas de Pinus cembroides Zucc.Rev. Chapingo Ser. Cienc. Forest. Ambiente12133–138.

  • Google Scholar

• 44

  GoodrichB. A.WaringK. M.AutyD.Sánchez MeadorA. J. (2018). Interactions of management and white pine blister rust on Pinus strobiformis regeneration abundance in southwestern USA.For. Int. J. For. Res.91492–505. 10.1093/forestry/cpy009

  • CrossRef
  • Google Scholar

• 45

  GoodrichB. A.WaringK. M.KolbT. E. (2016). Genetic variation in Pinus strobiformis growth and drought tolerance from southwestern US populations.Tree Physiol.361219–1235. 10.1093/treephys/tpw052

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  GregoriusH. R.BergmannF.WehenkelC. (2003). Analysis of biodiversity across levels of biological organization: a problem of defining traits.Perspect. Plant Ecol. Evol. Syst.5209–218. 10.1078/1433-8319-00035

  • CrossRef
  • Google Scholar

• 47

  GrimeJ. P.JeffreyD. W. (1965). Seedling establishment in vertical gradients of sunlight.J. Ecol.53621–642. 10.2307/2257624

  • CrossRef
  • Google Scholar

• 48

  Hernández-VelascoJ.Hernández-DíazJ. C.FladungM.Cañadas-LópezÁPrieto-RuízJ. ÁWehenkelC. (2017). Spatial genetic structure in four Pinus species in the Sierra Madre Occidental, Durango, México.Can. J. For. Res.4773–80. 10.1139/cjfr-2016-0154

  • CrossRef
  • Google Scholar

• 49

  HerreraC. M.CerdáX.GarcíaM. B.GuitiánJ.MedranoM.ReyP. J.et al (2002). Floral integration, phenotypic covariance structure and pollinator variation in bumblebee-pollinated Helleborus foetidus.J. Evol. Biol.15108–121. 10.1046/j.1420-9101.2002.00365.x

  • CrossRef
  • Google Scholar

• 50

  HerreraC. M.JordanoP.Lopez-SoriaL.AmatJ. A. (1994). Recruitment of a mast-fruiting, bird-dispersed tree: bridging frugivore activity and seedling establishment.Ecol. Monogr.64315–344. 10.2307/2937165

  • CrossRef
  • Google Scholar

• 51

  HiemstraP. H.PebesmaE. J.TwenhöfelC. J.HeuvelinkG. B. (2009). Real-time automatic interpolation of ambient gamma dose rates from the Dutch radioactivity monitoring network.Comput. Geosci.351711–1721. 10.1016/j.cageo.2008.10.011

  • CrossRef
  • Google Scholar

• 52

  HoffR.HagleS. (1990). “Diseases of whitebark pine with special emphasis on white pine blister rust,” in Proceedings of the Symposium on Whitebark Pine Ecosystems: Ecology and Management of a High-Mountain Resource, edsSchmidtW. C.McDonaldK. J. (Washington, DC: USDA), 179–190.

  • Google Scholar

• 53

  HutchinsonT. C. (1967). Comparative studies of the ability of species to withstand prolonged periods of darkness.J. Ecol.55291–299. 10.2307/2257878

  • CrossRef
  • Google Scholar

• 54

  IglesiasL. G.Solís-RamosL. Y.Viveros-ViverosH. (2012). Variación morfométrica en dos poblaciones naturales de Pinus hartwegii Lindl. del estado de Veracruz.Rev. Int. Bot. Exp.81239–246.

  • Google Scholar

• 55

  JarvisA.ReuterH. I.NelsonA.GuevaraE. (2008). Hole-Filled SRTM for the GLOBE VERSION 4, Available from the CGIAR-CSI SRTM 90m Database. Available online at: http://srtm.csi.cgiar.org(accessed August 26, 2019).

  • Google Scholar

• 56

  JiM.DengJ.YaoB.ChenR.FanZ.GuanJ.et al (2016). Ecogeographical variation of 12 morphological traits within Pinus tabulaeformis: the effects of environmental factors and demographic histories.J. Plant Ecol.10386–396. 10.1093/jpe/rtw033

  • CrossRef
  • Google Scholar

• 57

  JiM.ZhangX.WangZ.ZhangQ.DengJ. (2011). Intra-versus inter-population variation of cone and seed morphological traits of Pinus tabulaeformis Carr. in northern China: impact of climate-related conditions.Polish J. Ecol.59717–727.

  • Google Scholar

• 58

  JordanoP. (1995). Frugivore-mediated selection on fruit and seed size: birds and St. Lucie’s cherry, Prunus mahaleb.Ecology762627–2639. 10.2307/2265833

  • CrossRef
  • Google Scholar

• 59

  KaplanD. R. (2001). The science of plant morphology: definition, history, and role in modern biology.Am. J. Bot.881711–1741. 10.2307/3558347

  • CrossRef
  • Google Scholar

• 60

  KeaneR. E.CushmanS. A. (2018). Best Friends Forever: The Whitebark Pine and Clark’s Nutcracker. The Wildlife Professional.Washington, DC: USDA.

  • Google Scholar

• 61

  KralR. (1993). ““Pinus”,” in Flora of North America North of Mexico, ed.Flora of North America Editorial Committee (New York, NY: Oxford University Press), 373–398.

  • Google Scholar

• 62

  KruskalW. H.WallisW. A. (1952). Use of ranks in one-criterion variance analysis.J. Am. Stat. Assoc.47583–621. 10.1080/01621459.1952.10483441

  • CrossRef
  • Google Scholar

• 63

  KrzanowskiW. J. (1987). Selection of variables to preserve multivariate data structure, using principal components.Appl. Stat.3622–33. 10.2307/2347842

  • CrossRef
  • Google Scholar

• 64

  KuhnM. (2012). Variable Importance Using the Caret Package. Available online at: http://cran.cermin.lipi.go.id/web/packages/caret/vignettes/caretSelection.pdf(accessed September 25, 2019).

  • Google Scholar

• 65

  LandguthE. L.CushmanS. A. (2010). CDPOP: a spatially explicit cost distance population genetics program.Mol. Ecol. Resour.10156–161. 10.1111/j.1755-0998.2009.02719.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 66

  LandguthE. L.ForesterB. R.EckertA. J.ShirkA. J.MenonM.WhippleA.et al (2020). Modelling multilocus selection in an individual-based, spatially-explicit landscape genetics framework.Mol. Ecol. Resourc.20605–615. 10.1111/1755-0998.13121

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  LandguthE. L.HoldenZ. A.MahalovichM. F.CushmanS. A. (2017). Using landscape genetics simulations for planting blister rust resistant whitebark pine in the US Northern Rocky Mountains.Front. Genet.8:1–12. 10.3389/fgene.2017.00009

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 68

  Leal-SáenzA.WaringK. M.SniezkoR.MenonM.Hernández-DíazJ. C.López-SanchezC. A.et al (2020). Differences in cone and seed morphology of Pinus strobiformis and Pinus ayacahuite.Southwest. Nat. (in press).

  • Google Scholar

• 69

  LeishmanM. R.WestobyM. (1994b). The role of seed size in seedling establishment in dry soil conditions–experimental evidence from semi-arid species.J. Ecol.82249–258. 10.2307/2261293

  • CrossRef
  • Google Scholar

• 70

  LeishmanM. R.WestobyM. (1994a). The role of large seed size in shaded conditions: experimental evidence.Funct. Ecol.8205–214. 10.2307/2389903

  • CrossRef
  • Google Scholar

• 71

  LeishmanM. R.WrightI. J.MolesA. T.WestobyM. (2000). “Chapter 2: The evolutionary ecology of seed size,” in Seeds: the Ecology of Regeneration in Plant Communities, 2nd Edn, ed.FennerM. (Oxford: CABI Publishing), 31–57. 10.1079/9780851994321.0031

  • CrossRef
  • Google Scholar

• 72

  LeslieA. B.BeaulieuJ. M.MathewsS. (2017). Variation in seed size is structured by dispersal syndrome and cone morphology in conifers and other nonflowering seed plants.New Phytol.216429–437. 10.1111/nph.14456

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  LindB. M.MenonM.BolteC. E.FaskeT. M.EckertA. J. (2018). The genomics of local adaptation in trees: are we out of the woods yet?Tree Genet. Genomes14:29. 10.1007/s11295-017-1224-y

  • CrossRef
  • Google Scholar

• 74

  LooneyC. E.WaringK. M. (2013). Pinus strobiformis (southwestern white pine) stand dynamics, regeneration, and disturbance ecology: a review.For. Ecol. Manag.28790–102. 10.1016/j.foreco.2012.09.020

  • CrossRef
  • Google Scholar

• 75

  LordJ.EganJ.CliffordT.JuradoE.LeishmanM.WilliamsD.et al (1997). Larger seeds in tropical floras: consistent patterns independent of growth form and dispersal mode.J. Biogeogr.24205–211. 10.1046/j.1365-2699.1997.00126.x

  • CrossRef
  • Google Scholar

• 76

  MaronnaR. A.MartinR. D.YohaiV. J.Salibián-BarreraM. (2018). Robust Statistics: Theory and Methods (with R).New York, NY: Wiley.

  • Google Scholar

• 77

  MartínezI.GarcíaD.ObesoJ. R. (2007). Allometric allocation in fruit and seed packaging conditions the conflict among selective pressures on seed size.Evol. Ecol.21517–533. 10.1007/s10682-006-9132-x

  • CrossRef
  • Google Scholar

• 78

  MazerS. J. (1989). Ecological, taxonomic, and life history correlates of seed mass among indiana dune angiosperms: ecological archives M059-001.Ecol. Monogr.59153–175. 10.2307/2937284

  • CrossRef
  • Google Scholar

• 79

  McGarigalK.CushmanS.StaffordS. (2000). Multivariate Statistics for Wildlife and Ecology Research.New York: Springer Science.

  • Google Scholar

• 80

  McGarigalK.CushmanS. A. (2002). Comparative evaluation of experimental approaches to the study of habitat fragmentation effects.Ecol. Appl.12335–345. 10.1890/1051-07612002012[0335:CEOEAT]2.0.CO;2

  • CrossRef
  • Google Scholar

• 81

  MenonM.BagleyJ. C.FriedlineC. J.WhippleA. V.SchoettleA. W.Leal-SàenzA.et al (2018). The role of hybridization during ecological divergence of southwestern white pine (Pinus strobiformis) and limber pine (P. flexilis).Mol. Ecol.271245–1260. 10.1111/mec.14505

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 82

  MenonM.LandguthE.Leal-SaenzA.BagleyJ. C.SchoettleA. W.WehenkelC.et al (2020). Tracing the footprints of a moving hybrid zone under a demographic history of speciation with gene flow.Evol. Appl.13195–209. 10.1111/eva.12795

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 83

  Moreno-LetelierA.PiñeroD. (2009). Phylogeographic structure of Pinus strobiformis Engelm. across the Chihuahuan Desert filter-barrier.J. Biogeogr.36121–131. 10.1111/j.1365-2699.2008.02001.x

  • CrossRef
  • Google Scholar

• 84

  MorgensternE. K. (1996). Geographic Variation in Forest Trees. Genetics Basis and Application of Knowledge in Silviculture.Vancouver, BC: UBC Press.

  • Google Scholar

• 85

  NgF. S. P. (1978). “Strategies of establishment in Malayan forest trees,” in Tropical Trees as Living Systems, edsTomlinsonP. B.ZimmermanM. H. (Cambridge: Cambridge University Press), 129–162.

  • Google Scholar

• 86

  O’BrienR. A. (2002). Arizona’s Forest Resources, 1999. USDA Forest Service Resource Bulletin RMRS-RB-2. Ogden, UT: USDA.

  • Google Scholar

• 87

  OksanenJ.BlanchetF. G.KindtR.LegendreP.O’HaraR. B.SimpsonG. L.et al (2010). Vegan: Community Ecology Package. R Package Version 1.17-4.

  • Google Scholar

• 88

  PebesmaE. J.BivandR. S. (2005). Classes and methods for spatial data in R.R News59–13.

  • Google Scholar

• 89

  PlatenkampG. A.ShawR. G. (1993). Environmental and genetic maternal effects on seed characters in Nemophila menziesii.Evolution47540–555. 10.1111/j.1558-5646.1993.tb02112.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 90

  PriceG. R. (1970). Selection and covariance.Nature227520–521. 10.1038/227520a0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 91

  Prieto-RuízJ. A.Bustamante-GarcíaV.Muñoz-FloresH. J.Álvarez-ZagoyaR. (2014). “Capítulo 2. Análisis de conos y semillas en coníferas,” in Técnicas en el Manejo Sustentable de los Recursos Naturales. Cuerpo Académico “Manejo de Recursos Naturales y Sustentabilidad”, edsGarzaF.GuevaraO. J. A.VillalónH. M.CarrilloA. O. (Monterrey: Impresión, Universidad Autónoma de Nuevo León), 27–44.

  • Google Scholar

• 92

  R Development Core Team (2017). R: A Language and Environment for Statistical Computing.Vienna: R Foundation for Statistical Computing.

  • Google Scholar

• 93

  RawatK.BakshiM. (2011). Provenance variation in cone, seed and seedling characteristics in natural populations of Pinus wallichiana A.B. Jacks (Blue Pine) in India.Ann. For. Res.5439–55.

  • Google Scholar

• 94

  RoachD. A. (1987). Variation in seed and seedling size in Anthoxanthum odoratum.Am. Midland Natural.117258–264. 10.2307/2425967

  • CrossRef
  • Google Scholar

• 95

  SalasE. A. L.ValdezR.MichelS. (2017). Summer and winter habitat suitability of Marco Polo argali in southeastern Tajikistan: a modeling approach.Heliyon3:e00445. 10.1016/j.heliyon.2017.e00445

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 96

  SalisburyE. J. (1942). The Reproductive Capacity Of Plants: Studies in Quantitative Biology.London: G. Bell & Sons.

  • Google Scholar

• 97

  SalisburyE. J. (1974). Seed size and mass in relation to environment.Proc. R. Soc. Lond. Ser. B. Biol. Sci.18683–88. 10.1098/rspb.1974.0039

  • CrossRef
  • Google Scholar

• 98

  SamanoS.TombackD. F. (2003). Cone opening phenology, seed dispersal, and seed predation in southwestern white pine (Pinus strobiformis) in southern Colorado.Écoscience10319–326. 10.1080/11956860.2003.11682780

  • CrossRef
  • Google Scholar

• 99

  SavolainenO.PyhäjärviT.KnürrT. (2007). Gene flow and local adaptation in trees.Annu. Rev. Ecol. Evol. Syst.38595–619. 10.1146/annurev.ecolsys.38.091206.095646

  • CrossRef
  • Google Scholar

• 100

  SchimpfD. J. (1977). Seed weight of Amaranthus retroflexus in relation to moisture and length of growing season.Ecology58450–453. 10.2307/1935621

  • CrossRef
  • Google Scholar

• 101

  SchoettleA. W.JacobiW. R.WaringK. M.BurnsK. S. (2018). Regeneration for resilience framework to support regeneration decisions for species with populations at risk of extirpation by white pine blister rust.New For.5089–114. 10.1007/s11056-108-9679-8

  • CrossRef
  • Google Scholar

• 102

  ShimadaT.TakahashiA.ShibataM.YagihashiT. (2015). Effects of within-plant variability in seed weight and tannin content on foraging behaviour of seed consumers.Funct. Ecol.291513–1521. 10.1111/1365-2435.12464

  • CrossRef
  • Google Scholar

• 103

  ShirkA. J.CushmanS. A.WaringK. M.WehenkelC.Leal-SáenzA.ToneyC.et al (2018). Southwestern white pine (Pinus strobiformis) species distribution models project a large range shift and contraction due to regional climatic changes.For. Ecol. Manag.411176–186. 10.1016/j.foreco.2018.01.025

  • CrossRef
  • Google Scholar

• 104

  SniezkoR. A.KegleyA. J.DanchokR. (2008). White pine blister rust resistance in North American, Asian and European species - results from artificial inoculation trials in Oregon.Ann. For. Res.5153–66.

  • Google Scholar

• 105

  SorensenF. C.MilesR. S. (1978). Cone and seed weight relationships in Douglas-fir from Western and Central Oregon.Ecology59641–644. 10.2307/1938763

  • CrossRef
  • Google Scholar

• 106

  SorkV. L.AitkenS. N.DyerR. J.EckertA. J.LegendreP.NealeD. B. (2013). Putting the landscape into the genomics of trees: approaches for understanding local adaptation and population responses to changing climate.Tree Genet. Genomes9901–911. 10.1007/s11295-013-0596-x

  • CrossRef
  • Google Scholar

• 107

  SteinhoffR. J.AndresenJ. W. (1971). Geographic variation in Pinus flexilis and Pinus strobiformis and its bearing on their taxonomic status.Silvae Genet.20159–167.

  • Google Scholar

• 108

  StrombergJ. C.PattenD. T. (1990). Variation in seed size of a southwestern riparian tree, Arizona Walnut (Juglans major).Am. Midland Natural.124269–277. 10.2307/2426176

  • CrossRef
  • Google Scholar

• 109

  TeleniusA.TorstenssonP. (1991). Seed wings in relation to seed size in the genus Spergularia.Oikos61216–222. 10.2307/3545339

  • CrossRef
  • Google Scholar

• 110

  TombackD. F.AchuffP. (2010). Blister rust and western forest biodiversity: ecology, values and outlook for white pines.For. Pathol.40186–225. 10.1111/j.1439-0329.2010.00655.x

  • CrossRef
  • Google Scholar

• 111

  Van Den WollenbergA. L. (1977). Redundancy analysis an alternative for canonical correlation analysis.Psychometrika42207–219. 10.1007/BF02294050

  • CrossRef
  • Google Scholar

• 112

  VenablesW. N.RipleyB. D. (1999). “Chapter 10: Tree-based methods,” in Modern Applied Statistics with S-PLUS, 3rd Edn, edsChambersJ.EddyW.HärdleW.SheatherS.TierneyL. (New York, NY: Springer-Verlag Press), 303–327.

  • Google Scholar

• 113

  Villagómez-LozaM. A.Bello-GonzálezM. A.Isarain-ChávezE. (2014). Pinus strobiformis Engelmann: Nueva localidad para Guanajuato, México.Agrociencia48615–625.

  • Google Scholar

• 114

  Villalobos-ArámbulaA. R.Pérez de la RosaJ. P.AriasA.RajoraO. P. (2014). Cross-species transferability of eastern white pine (Pinus strobus) nuclear microsatellite markers to five Mexican white pines.Genet. Mol. Res.137571–7576. 10.4238/2014.september.12.24

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 115

  ViolleC.NavasM. L.VileD.KazakouE.FortunelC.HummelI.et al (2007). Let the concept of trait be functional!Oikos116882–892. 10.1111/j.2007.0030-1299.15559.x

  • CrossRef
  • Google Scholar

• 116

  VittinghoffE.McCullochC. E. (2007). Relaxing the rule of ten events per variable in logistic and Cox regression.Am. J. Epidemiol.165710–718. 10.1093/aje/kwk052

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 117

  WahidN.González-MartínezS. C.El HadramiI.BoulliA. (2006). Variation of morphological traits in natural populations of maritime pine (Pinus pinaster Ait.) in Morocco.Ann. For. Sci.6383–92. 10.1051/forest:20050100

  • CrossRef
  • Google Scholar

• 118

  WangB.IvesA. R. (2017). Tree-to-tree variation in seed size and its consequences for seed dispersal versus predation by rodents.Oecologia183751–762. 10.1007/s00442-016-3793-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 119

  WehenkelC.Mariscal-LuceroS. R.Jaramillo-CorreaJ. P.López-SánchezC. A.Vargas-HernándezJ. J.Sáenz-RomeroC. (2017). “Genetic diversity and conservation of Mexican forest trees,” in Biodiversity and Conservation of Woody Plants. Sustainable Development and Biodiversity, Vol. 17edsAhujaM.JainS. (Cham: Springer), 37–67. 10.1007/978-3-319-66426-2_2

  • CrossRef
  • Google Scholar

• 120

  WehenkelC.Quiñones-PérezC. Z.Hernández-DíazJ. C.López-SánchezC. A. (2014). “Ecology and Genetics of Pinus strobiformis,” in Proceedings of the IUFRO Joint Conference: Genetics of Five-Needle Pines, Rusts of Forest Trees, and Strobusphere, edsSchoettleA. W.SniezkoR. A. (Fort Collins, CO: US Department of Agriculture, Forest Service, Rocky Mountain Research Station), 1–6.

  • Google Scholar

• 121

  WheelerN. C.GuriesR. P. (1982). Population structure, genic diversity, and morphological variation in Pinus contorta Dougl.Can. J. For. Res.12595–606. 10.1139/x82-091

  • CrossRef
  • Google Scholar

• 122

  WickhamH.ChangW.WickhamM. H. (2013). Package ‘ggplot2’. Computer Software Manual]. Available online at: http://cran.r-project.org/web/packages/ggplot2/ggplot2.pdf(accessed September 4, 2019).

  • Google Scholar

• 123

  WilliamsC. K.EngelhardtA.CooperT.MayerZ.ZiemA.ScruccaL.et al (2018). Package ‘caret’, version 6.0-80. [Online]. Available online at: https://cran.r-project.org/web/packages/caret/caret.pdf(accessed September 26, 2019).

  • Google Scholar

• 124

  WinnA. A. (1988). Ecological and evolutionary consequences of seed size in Prunella vulgaris.Ecology691537–1544. 10.2307/1941651

  • CrossRef
  • Google Scholar

• 125

  ZhangC.LiX.ChenL.XieG.LiuC.PeiS. (2016). Effects of topographical and edaphic factors on tree community structure and diversity of subtropical mountain forests in the lower lancang river basin.Forests7:222. 10.3390/f7100222

  • CrossRef
  • Google Scholar