Intraspecific Niche Models for Ponderosa Pine (Pinus ponderosa) Suggest Potential Variability in Population-Level Response to Climate Change

Kaitlin C Maguire; Douglas J Shinneman; Kevin M Potter; Valerie D Hipkins

doi:10.1093/sysbio/syy017

Intraspecific Niche Models for Ponderosa Pine (Pinus ponderosa) Suggest Potential Variability in Population-Level Response to Climate Change

Maguire, Kaitlin C; Shinneman, Douglas J; Potter, Kevin M; Hipkins, Valerie D 2018-11-01 00:00:00 Abstract Unique responses to climate change can occur across intraspecific levels, resulting in individualistic adaptation or movement patterns among populations within a given species. Thus, the need to model potential responses among genetically distinct populations within a species is increasingly recognized. However, predictive models of future distributions are regularly fit at the species level, often because intraspecific variation is unknown or is identified only within limited sample locations. In this study, we considered the role of intraspecific variation to shape the geographic distribution of ponderosa pine (Pinus ponderosa), an ecologically and economically important tree species in North America. Morphological and genetic variation across the distribution of ponderosa pine suggest the need to model intraspecific populations: the two varieties (var. ponderosa and var. scopulorum) and several haplotype groups within each variety have been shown to occupy unique climatic niches, suggesting populations have distinct evolutionary lineages adapted to different environmental conditions. We utilized a recently available, geographically widespread dataset of intraspecific variation (haplotypes) for ponderosa pine and a recently devised lineage distance modeling approach to derive additional, likely intraspecific occurrence locations. We confirmed the relative uniqueness of each haplotype-climate relationship using a niche-overlap analysis, and developed ecological niche models (ENMs) to project the distribution for two varieties and eight haplotypes under future climate forecasts. Future projections of haplotype niche distributions generally revealed greater potential range loss than predicted for the varieties. This difference may reflect intraspecific responses of distinct evolutionary lineages. However, directional trends are generally consistent across intraspecific levels, and include a loss of distributional area and an upward shift in elevation. Our results demonstrate the utility in modeling intraspecific response to changing climate and they inform management and conservation strategies, by identifying haplotypes and geographic areas that may be most at risk, or most secure, under projected climate change. Ecological niche modeling, Pinus ponderosa, intraspecific variation, haplotype, range shift, conservation, climate-plant relationships A species’ response to environmental change is not always uniform across its distribution as genotypic and phenotypic differences within a species can interact with local environmental conditions and result in local adaptation (Valladares et al. 2014). Even though most biogeographical modeling studies treat a species as a single unit, its response to climate change will typically vary across space, either through plastic gene expression, adaptation, extirpation, or migration (Sork et al. 2010; Oney et al. 2013; Valladares et al. 2014). Therefore, population and phylogeographic structure should not be ignored when modeling species distributions, inferring niche structure, or predicting how climate will affect species distributions (Pfenninger et al. 2007; Pearman et al. 2010). Not including phylogeographic structure may lead to underrepresentation of certain lineages in species models (Pearman et al. 2010), especially when they are clearly defined and occupy distinct niche spaces (Hällfors et al. 2016). Distinct lineages defined by genetic variation are the units of evolutionary change and thus crucial for providing insight into the ability of a species to adapt to environmental change. Predicting intraspecific geographic response to environmental factors is becoming increasingly feasible as research reveals genetic variation across the distribution of individual species (e.g., Pearman et al. 2010; Espíndola et al. 2012; D’Amen et al. 2013; Yannic et al. 2014; Hällfors et al. 2016; Prasad and Potter 2017). Modeling the ecological niches of populations or lineages may be especially important for species with spatial variation in genotypic or phenotypic traits that suggest they are in the process of differentiating into two or more species (Valladares et al. 2014). Ponderosa pine (Pinus ponderosa) is a widespread, well-studied, and economically important species in Western North America, and its substantial morphological variation (Wells 1964) suggests it may be in early stages of differentiation into multiple species (Wang 1977, Jaramillo-Correa et al. 2009). Two varieties within the species are recognized—the Pacific variety, P. ponderosa var. ponderosa Laws., has open, plume-like foliage, and a low proportion of two-needle fascicles as opposed to the Rocky Mountain variety, P. ponderosa var. scopulorum Engelm., which has compact and brush-like foliage and moderate to high proportions of two-needle fascicles. The separation of the two varieties is also established by leaf oil terpene composition (Rudloff and Lapp 1992), mitochondrial DNA (Johansen and Latta 2003; Potter et al. 2013), highly polymorphic nuclear microsatellite and isozyme markers (Potter et al. 2015), and allozymes (Niebling and Conkle 1990). The two varieties have similar temperature limitations but occur within different precipitation regimes, with var. ponderosa locations generally corresponding to higher winter moisture and var. scopulorum to higher summer moisture availability (Norris et al. 2006; Shinneman et al. 2016). It is thought that the varieties diverged at least before the last glacial maximum, and more likely around 250,000 years ago (Lascoux et al. 2004), with Pleistocene glacial and interglacial cycles leading to further genetic diversification through isolation in localized refugia, followed by subsequent migration, hybridization, and introgression (Shinneman et al. 2016). Within each variety, several geographically structured haplotypes have also been identified using mitochondrial DNA (Johansen and Latta 2003; Potter et al. 2013) and plastome sequences (Wofford et al. 2014), and genetic clusters have been identified from nuclear DNA microsatellite markers (Potter et al. 2015). Geographic structure identified from mtDNA suggests each haplotype resulted from unique patterns of migration, isolation, and local adaptation (Potter et al. 2013). Mitochondrial DNA is maternally inherited and dispersed via seeds with limited dispersal compared with pollen-dispersed chloroplast DNA, and geographic structure in genetic differentiation is retained longer (Neale and Sederoff 1989; Petit et al. 1993). Thus, the haplotypes may represent evolutionarily distinct units capable of responding differently to climate change as a result of unique adaptations to refugial climate conditions and long-term isolation during glacial periods (Potter et al. 2013, Shinneman et al. 2016). Given unique intraspecific lineages and the potential for distinctive climate niches, the need to model distributions of ponderosa pine at the subspecific level has been recognized (Norris et al. 2006; Rehfeldt et al. 2014, Shinneman et al. 2016). Shinneman et al. (2016) modeled the contemporary and historical climate niche distributions of the two varieties, as well as 10 geographically structured mtDNA haplotypes identified by Potter et al. (2013) to further infer the evolutionary and phylogeographic history of these lineages. Recent modeling research has also suggested that ponderosa pine may be threatened by changing climate, but that the response may differ for the two varieties (Rehfeldt et al. 2014). However, the potential for future climate change to shape the distribution of ponderosa pine populations, as defined by haplotypic diversity, has not been explored, and has only been studied for a handful of other species relative to intraspecific variation (Pearman et al. 2010; Bálint et al. 2011; Bentio Garzón et al. 2011; Oney et al. 2013; Prasad and Potter 2017). Given past climate influences on ponderosa pine evolutionary trajectories and genetic lineages, we sought to better understand potential future distributions of the species by specifically considering responses to expected climate change of the two well-known and broadly distributed varieties and the haplotypes associated with them. Specifically, we were interested in determining how haplotype-level ecological niche models (ENMs) might differ from variety-level ENMs in terms of predicting recent and future geographic ranges in relation to climate, and exploring the advantages and challenges of estimating ENMs below the traditional level of species. To accomplish this objective, we employed a three-tiered approach. First, to address potentially inadequate geographic coverage of samples across populations, we utilized a lineage distance modeling approach (Rosauer et al. 2015) to derive likely haplotype occurrence locations (presences) by leveraging distribution information for the two varieties and known haplotype locations. Second, both actual and derived presence-absence datasets were employed in an ensemble of ENMs to project recent (1961–1990) and future (2060, 2090) distributions for ponderosa pine across three intraspecific levels: 1) for each individual haplotype; 2) for all haplotypes within a variety combined; and 3) for each corresponding variety. Third, to quantify the uniqueness of each estimated haplotype niche space relative to all others, we applied a post hoc, niche-overlap analysis. The combination of these methodological approaches allowed us to provide relatively robust projections of future distributions of ponderosa pine across intraspecific levels, and permitted a more nuanced interpretation of potentially unique population-level responses to climate change that may better inform conservation efforts for the species. Methods Haplotype Occurrence Data We examined eight of the ten haplotypes identified in Potter et al. (2013). We did not examine haplotypes 9 and 10, each limited to a small area of California, because of restricted and localized populations less amenable to niche modeling. Haplotypes were identified from the mtDNA nad1 second intron minisatellite region of foliage samples collected from 3113 trees representing 104 populations of ponderosa pine across its range (Potter et al. 2013). At least 20 individual trees were sampled in each population (except two), with most populations encompassing at least 30 sampled trees. Sampled trees were at least 100 m apart to encompass the entire range of the population’s genetic composition. Haplotypes 1, 5, and 8 are within the distribution of var. ponderosa and group together in phylograms of genetic distance (Potter et al. 2013). Haplotypes 2, 3, 4, 6, and 7 are within the distribution of var. scopulorum and are phylogenetically distinct from the haplotypes in var. ponderosa (Potter et al. 2013). Duplicate occurrences of the same haplotype that fell within the same geographic grid cell (800 m |$\times$| 800 m) were removed from the dataset. Climate Data Climate variables derived using thin-plate splines to interpolate 1961–1990 monthly temperature and precipitation data at 800 m resolution were obtained from the Research on Forest Climate Change website (charcoal.cnre.vt.edu/climate; Rehfeldt 2006; Crookston and Rehfeldt 2008), as well as an ensemble of 17 global climate models (Supplementary Appendix S1 available on Dryad at http://dx.doi.org/10.5061/dryad.jj264) for the RCP60 (medium-high greenhouse gas concentration pathway) scenario for 2060 and 2090. These climate layers were restricted to the western United States. Additional climate variables were derived from these climate variables (Table 1). We removed one of any two variables that was highly correlated (Pearson correlation |$>$| 0.7), keeping the more biologically relevant variable of a correlated pair (Rehfeldt et al. 2014; Shinneman et al. 2016). This resulted in seven variables used to fit the final models (except for random forest models of the varieties—see below): growing season precipitation, summer-spring precipitation balance, degree days below 0|$^{\circ}$|C, mean temperature of the warmest month, summer-winter temperature differential, annual moisture index, and precipitation ratio. Variables were log or square root transformed if needed to normalize data. Table 1. Climate variables considered for niche modeling of ponderosa pine varieties and haplotypes Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| MAP |$=$| mean annual precipitation; GSP |$=$| growing season precipitation; MTCM |$=$| mean temperature of the coldest month; MTWM |$=$| mean temperature of the warmest month. |$^{\mathrm{a}}$|All variables were used in an initial random forest model for variety distributions, and seven variables with Pearson correlations |$<$|0.7 were used in final niche models. |$^{\mathrm{b}}$|Derived variables are indicated. Open in new tab Table 1. Climate variables considered for niche modeling of ponderosa pine varieties and haplotypes Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| MAP |$=$| mean annual precipitation; GSP |$=$| growing season precipitation; MTCM |$=$| mean temperature of the coldest month; MTWM |$=$| mean temperature of the warmest month. |$^{\mathrm{a}}$|All variables were used in an initial random forest model for variety distributions, and seven variables with Pearson correlations |$<$|0.7 were used in final niche models. |$^{\mathrm{b}}$|Derived variables are indicated. Open in new tab Lineage Distance Modeling for Haplotypes Although predicting the niches of intraspecific populations is more feasible given the rise in research that uncovers genetic variation of species across their distributions, often there are not enough sampled occurrences of each population to define their climate niche space and fit robust ENMs, limiting the ability to hindcast or forecast their distributions (Rosauer et al. 2015). In addition, it may not be appropriate to fit an ENM directly to only known haplotype occurrences because these occurrences may be governed by other factors than environmental conditions, such as competition. In the case of ponderosa pine, there may be too few sampled occurrences among intraspecific populations to produce robust ENMs, and sampled occurrences alone most likely do not represent the entire climatic niche of each haplotype. Thus, to overcome low sample sizes of known haplotype occurrences, we employed a recently developed method, lineage distance modeling (Rosauer et al. 2015), that probabilistically estimates the likely distribution of intraspecific lineages, therefore, increasing the number of occurrence points to fit ENMs. This is accomplished by first fitting an initial ENM for a species, and then assigning pixels within the species distribution to haplotypes based on the distance to known haplotype occurrences, taking into account barriers to dispersal using a measured cost distance (Fig. 1). The cost-distance analysis is calculated as |$-$|log(habitat suitability) where habitat suitability is the probability of occurrence from the species ENM. The output is a probability map for each lineage, such that all estimated probabilities for a cell sum to the total probability of the species occupying the cell from the original ENM. For this analysis, ENMs were fit for the two varieties separately (as described below), rather than the entire species, because they have distinct climate niches (Norris et al. 2006). The resulting distribution maps for each haplotype are effectively defined by the niche of the variety and geographically constrained by theoretical dispersal limitations. Custom python code from Rosauer et al. (2015) was used in ArcMap 10.4.1 to develop the haplotype probability maps. Estimated probability maps for each haplotype were compared with known haplotype locations using the Brier Score (Brier 1950), a measure of the accuracy of probabilistic predictions to binary (presence/absence) observations calculated as the mean squared difference between the two. The lower the Brier Score, the better the predictions. Figure 1. Open in new tabDownload slide Workflow of lineage distance modeling and ENM shown for var. ponderosa and haplotypes 1, 5, and 8. First, an ENM is fit for the variety using RF (Step 1). Then, lineage distance modeling is used to produce distribution maps for each haplotype (Step 2). These distribution maps are sampled, and ENMs are fit for each haplotype using RF, MARS, and NPMR (Step 3). Lastly, ENMs are fit for the variety using RF, MARS, and NPMR (Step 4). An ensemble distribution is created for each haplotype and the variety by averaging the RF, MARS, and NPMR models. Figure 1. Open in new tabDownload slide Workflow of lineage distance modeling and ENM shown for var. ponderosa and haplotypes 1, 5, and 8. First, an ENM is fit for the variety using RF (Step 1). Then, lineage distance modeling is used to produce distribution maps for each haplotype (Step 2). These distribution maps are sampled, and ENMs are fit for each haplotype using RF, MARS, and NPMR (Step 3). Lastly, ENMs are fit for the variety using RF, MARS, and NPMR (Step 4). An ensemble distribution is created for each haplotype and the variety by averaging the RF, MARS, and NPMR models. Variety and Haplotype Niche Modeling To fit ENMs for each variety, 20,000 cells were subsampled from ponderosa pine distribution maps of the USDA Forest Service’s National Individual Tree Species Atlas (https://www.fs.fed.us/foresthealth/applied_sciences/mapping-reporting/remote-sensing/index.shtml) that were also occupied by vegetation types that contain ponderosa pine in the U.S. Geological Survey National Gap Analysis Program (https://gapanalysis.usgs.gov) map (Fig. 1). Because both of these mapped sources were derived from classifications of satellite imagery, we used agreement between the two datasets for presence locations to minimize errors of commission (similar to Shinneman et al. 2016). Each of the occupied cells was assigned to one of the two varieties based on its geographic position. In transition zones between the two varieties (southern California and central Montana), the closest known haplotype occurrence was used to assign cells to a variety. Lastly, known haplotype occurrences were added to their corresponding variety. In the case of multiple occurrences of a single haplotype at a given location, only one record per that haplotype per grid cell (800 m |$\times$| 800 m) was kept for analysis. An initial 20,000 absence points were randomly subsampled, and 500 additional absences were generated for each variety from the cells occupied by the other variety. To fit ENMs for each haplotype, original haplotype presences from Potter et al. (2013) were augmented with up to 1000 cells sampled from the distribution map derived from the lineage distance modeling method (Fig. 1). Sampling from the lineage distance map was weighted so that 75% of samples had probabilities between 0.9 and 1, and 25% were between 0.5 and 0.89. Although this weighting scheme is a conservative approach to estimating haplotype occurrences and minimizes potential overestimation in peripheral populations, the lineage distance modeling also generally added samples in peripheral locations relative to the actual haplotype locations (Supplementary Appendix S2 available on Dryad), potentially balancing out any bias against peripheral populations. Indeed, after testing two less conservative weighting approaches, we found no difference between the ability of each sampling scheme to predict known haplotype locations (i.e., Brier Scores were not significantly different). Moreover, we also explored the possibility that using raw probability of occurrence values from the lineage distance modeling (analogous to using abundance) would produce less-restrictive ENMs for haplotypes, but the projected distributions were similar (or even more restricted, see Supplementary Appendix S3 available on Dryad). We also found little difference in Brier Scores using raw probability of occurrence values compared with the threshold sampling procedure (Supplementary Appendix S3 available on Dryad). Absences for each haplotype consisted of original absences from Shinneman et al. (2016) and absences generated in the previous step for the variety to which it belongs. We used the absences of each variety rather than of the entire species because the delineation of haplotypes between the varieties is robust based on mtDNA and nuclear DNA analyses (Potter et al. 2013; Potter et al. 2015). For all haplotype and variety niche models, data were split into 10 random partitions of training (80%) and testing (20%) sets. To minimize uncertainty in model algorithm choice, three different ENM modeling approaches were used to fit each partition (i.e., resulting in 30 models per haplotype/variety): random forest (RF), multivariate adaptive regression spline (MARS), and non-parametric multiplicative regression (NPMR). RF and MARS models were fit and projected in R v 3.3.1 (R Core Team 2014) using the randomForest (Liaw and Wiener 2002) and earth (Milborrow 2016) packages, respectively. NPMR models were run in HyperNiche 2.28 (McCune and Mefford 2009). All 23-climate variables were input into the initial RF models built for the varieties that provided inputs for the lineage distance modeling to explore all climate variables and minimize under-prediction of variety distributions (RF randomly chooses predictor variables and eliminates those that are not improving model performance; therefore, removal of correlated variables a priori is not needed.) However, to facilitate compatibility of results among different modeling approaches, final variety and haplotype RF models were fit with the same seven uncorrelated climate variables used in the MARS and NPMR models (Table 1). Each individual haplotype and variety model was projected onto current and future climatic scenarios, and an ensemble model was created by first averaging the 10 model runs of ecological niche maps estimated from each modeling approach, and then averaging the three resulting maps per model approach for each climate scenario. We also combined the distributions of haplotypes within each variety (H1, H5, and H8 for var. ponderosa; H2, H3, H4, H6, and H7 for var. scopulorum) at each time period in order to observe how the distribution of the variety produced from ENM haplotypes compared with the distribution produced by an ENM fit with variety occurrences. These combined distribution maps represent the probability of at least one of the haplotypes occurring in a particular cell (i.e., by calculating |$1 - X$|⁠, where |$X$| is the probability that all haplotypes are absent). Model performance for each ENM (haplotypes and varieties) was tested using the area under the receiver operating curve (AUC), a metric used when comparing continuous outputs (ENM occurrence probabilities) to binary observations (presence/absence dataset sampled from the lineage distance method), ranging from 0 to 1 in which a value of 1 indicates perfect agreement and a value of 0.5 indicating model performance is no better than random. We also evaluated model performance by comparing the known haplotype locations to the predicted occurrence probabilities using the Brier Score. Niche-Overlap Analyses The amount of niche overlap and niche similarity between haplotypes was calculated following Broennimann et al. (2012), using the “PCA-env option”, in which a Kernel smoother is applied to densities of species occurrence in a gridded environmental space calibrated on the available environmental space. The D metric (Schoener 1970) measures the amount of niche overlap between two haplotypes in the gridded environmental space (D |$=$| 1: complete overlap; D |$=$| 0: no overlap): \[ D=1-0.5 \left(\sum_{ij} \left| {z1}_{ij}-{z2}_{ij} \right|\right) \] where |${z1}_{ij}$| and |${z2}_{ij}$| are the occupancies of entity 1 and 2, respectively, and |$i$| and |$j$| refer to the cell corresponding to the |$i$|th and |$j$|th bins of the environmental variables. The similarity test measures the amount of overlap of the niche of one haplotype with a random subset in the second haplotype’s range, and vice versa, for 100 iterations, and compares this distribution of overlap to the observed D metric. The similarity test can indicate if two niches are significantly similar, significantly dissimilar, or nonsignificant (indicating not enough statistical power to determine either/or). Niche space here was defined using the same presence points used in ENM modeling. Niche overlap and similarity analyses were performed using the ecospat package (Broennimann et al. 2016) in R. Results Lineage Distance Modeling for Haplotypes Distribution maps of each haplotype from the lineage distance modeling are generally in agreement with known locations of each haplotype (Brier Scores |$<$| 0.2; Table 2; Supplementary Appendix S2 available on Dryad). However, there were limitations to this approach for certain haplotypes. Specifically, because the distributions of H2, H4, and H7 are geographically restricted (Fig. 3 in Supplementary Appendix S2 available on Dryad), fewer than 1000 occurrence points from the lineage distance maps for these haplotypes were available for ENM construction (Table 2; H4: |$n =$| 102, H2: |$n =$| 296, and H7: |$n =$| 253). Furthermore, because the ENM for var. scopulorum predicted only low probabilities of occurrence (⁠|$<$| 0.50) for the small and isolated populations of H2 in southern California/Nevada and H7 in Nevada, in those specific areas only the original haplotype locations (from Potter et al. 2013) were used to construct ENMs. Figure 3. Open in new tabDownload slide Percent area change and mean elevation change from the recent period to 2060 and 2090 for the varieties (Var.), combined haplotypes (Haps.), and individual haplotypes. Change values are calculated among time periods using cells with probability of occurrence estimates |$>$| 0.50. Figure 3. Open in new tabDownload slide Percent area change and mean elevation change from the recent period to 2060 and 2090 for the varieties (Var.), combined haplotypes (Haps.), and individual haplotypes. Change values are calculated among time periods using cells with probability of occurrence estimates |$>$| 0.50. Table 2. Lineage distance modeling and ensemble ENM evaluation Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Notes: Brier scores represent the mean squared difference between the predicted probabilities of occurrence and the observed binary presence/absence of each haplotype for the lineage distance method and ensemble ENM. Mean AUC values across 10 model runs for all three models for each variety and haplotype indicate how well ENMs discriminates between presences and absences. Lastly, the number of presences for each variety and haplotype used to build the ENMs are listed, with the number of actual locations (from Potter et al. 2013) followed by estimated locations (using the linear distance method) in parentheses. Open in new tab Table 2. Lineage distance modeling and ensemble ENM evaluation Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Notes: Brier scores represent the mean squared difference between the predicted probabilities of occurrence and the observed binary presence/absence of each haplotype for the lineage distance method and ensemble ENM. Mean AUC values across 10 model runs for all three models for each variety and haplotype indicate how well ENMs discriminates between presences and absences. Lastly, the number of presences for each variety and haplotype used to build the ENMs are listed, with the number of actual locations (from Potter et al. 2013) followed by estimated locations (using the linear distance method) in parentheses. Open in new tab Variety and Haplotype Niche Models ENMs for varieties and haplotypes were able to discriminate well between occupied and unoccupied cells (AUC scores |$>$| 0.9; Table 2), and predicted the occurrence of known haplotype occurrences well (Brier Score |$< = 2.0$|⁠, Table 2). Variable importance differed for each variety and haplotype as well as among model type (Supplementary Appendix S4 available on Dryad). Var. ponderosa.— The ensemble ENM for the variety matches its known distribution, although it slightly over-predicts in portions of the Great Basin (in Nevada, Southern Idaho, and Eastern Oregon) and the Wasatch Range in Northern Utah (Fig. 2a). Projections of the ENM for the variety to 2060 and 2090 (Fig. 2b and c) suggest a range contraction will occur, with a projected 24% loss in total area within high occurrence probability areas (⁠|$>$| 0.5 occurrence probability) by 2090 (compared with the recent period) (Fig. 3), especially at lower elevations in central Oregon, northeastern Washington, and northeastern California. This trend is compensated slightly with potential range extension into higher elevations, especially in central Idaho and the southern Sierra Nevada of California. Indeed the mean elevation for cells with high probability (⁠|$>$| 0.5) of occurrence is projected to increase by 132 m on average by 2090 (compared with the recent period) (Fig. 3). Moreover, the mean probability of all occupied cells with |$>$| 0.5 probability in the recent period is generally projected to decrease by 0.02 and then by 0.53 in 2060 and 2090, respectively (Fig. 4). The combined distribution for the three haplotypes (H1, H5, and H8) within the variety is not as extensive compared with the variety model for the recent climate period (Fig. 2d), but similar distributional trends, such as loss in central Oregon, are projected through 2090 (Fig. 2e–f). However, for the combined haplotype distribution there is no loss in total area and only a 57 m increase in mean elevation for cells with high probability (⁠|$>$| 0.5) of occurrence (Fig. 3), although mean probability values for all occupied cells in the recent period with |$>$| 0.5 probability occurrence are projected to decline for both 2060 and 2090, by 0.17 and 0.53, respectively (Fig. 4). Individual haplotype models generally project similar trends to the variety (Supplementary Appendix S5 available on Dryad), as highlighted here for H1 (Fig. 2g–i), which is projected to have a decrease in area (41%), an increase in mean elevation (132 m), and a decrease in mean probability (by 0.63, of occupied cells in the recent period with |$>$| 0.5 probability) for 2090 (Figs. 3 and 4). Figure 2. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. ponderosa (a–c), combined haplotype (H1, H5, and H8) models (d–f), and a single haplotype (H1) model (g–i). Predicted distributions in all cases are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. ponderosa (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Figure 2. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. ponderosa (a–c), combined haplotype (H1, H5, and H8) models (d–f), and a single haplotype (H1) model (g–i). Predicted distributions in all cases are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. ponderosa (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Figure 4. Open in new tabDownload slide Mean occurrence probability of occupied cells (with |$>$| 0.5 occurrence probability) from ensemble model of recent distributions of ponderosa pine lineages, combined haplotypes, and individual haplotypes and the mean occurrence probability of those cells in 2060 and 2090. Figure 4. Open in new tabDownload slide Mean occurrence probability of occupied cells (with |$>$| 0.5 occurrence probability) from ensemble model of recent distributions of ponderosa pine lineages, combined haplotypes, and individual haplotypes and the mean occurrence probability of those cells in 2060 and 2090. Var. scopulorum.— The modeled niche distribution of var. scopulorum matches the known distribution of the variety, although there is some over-prediction of the distribution in Montana and under prediction in the Great Plains (Fig. 5a). In general, the ensemble ENM of the variety and the individual haplotypes (H3 and H6) could not accurately model the northern portion of the variety’s distribution. The easternmost populations of the variety along the Niobrara River in Nebraska are not predicted by the ensemble ENMs most likely because these populations exist in microclimates along the river that are not available in surrounding areas, making it difficult for coarse scale climate data (800 m x 800 m resolution) to capture their presence. Figure 5. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. scopulorum (a–c), combined haplotype (H2, H3, H4, H6, and H7) models (d–f), and a single haplotype (H3) model (g–i). Predicted niche distributions are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. scopulorum (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Figure 5. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. scopulorum (a–c), combined haplotype (H2, H3, H4, H6, and H7) models (d–f), and a single haplotype (H3) model (g–i). Predicted niche distributions are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. scopulorum (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Similar to var. ponderosa, ENM projections for var. scopulorum in the future suggest a range contraction, with a 59% area loss and a 670 m increase in mean elevation by 2090 compared with the recent time period for areas with |$>$| 0.5 occurrence probability (Fig. 3). There is some predicted expansion northward into central Utah and western Colorado, and some expansion into southwest Montana even though there is predicted loss for most Montana populations. Occupied cells with probabilities |$>$| 0.5 in the recent period decline by 0.17 in 2060 and by 0.57 in 2090 (Fig. 4). The combined predicted niche distributions of the individual haplotypes within the variety (H2, H3, H4, H6, and H7) are more delimited than the variety (similar to the pattern observed with var. ponderosa haplotypes), but they have a similar pattern of loss and upward range shift as the variety, with a 23% loss in area and a mean elevational increase of 363 m for high probability (⁠|$>$| 0.5) occurrences by 2090 (Fig. 3). The probability of the combined haplotypes occurring in predicted occupied sites in the recent time period with |$>$| 0.5 probability in the recent time period decreases by 0.61 in 2090 (Fig. 4). Individual haplotypes within var. scopulorum follow a similar trend (Supplementary Appendix S5 available on Dryad), as highlighted by H3 (Fig. 5g–i), which is projected to have a decrease in area (40%), an increase in mean elevation (422 m), and a decrease in mean probability (by 0.63, of occupied cells in the recent period with |$>$| 0.5 probability) for 2090 (Figs. 3 and 4). Niche Overlap Haplotypes of ponderosa pine occupy distinct niche spaces in western North America, as no two haplotypes had statistically similar niche spaces (Similarity Test |$P \,{<}\,0.01$|⁠; Supplementary Appendix S6 available on Dryad). Thus, even though the niche spaces of some haplotypes overlap (e.g., H3 and H7), they are not significantly similar (Fig. 6 and Supplementary Appendix S6 available on Dryad). Haplotypes with larger distributional ranges tend to occupy more niche space (e.g., H1 vs. H4; Fig. 6; Supplementary Appendix S6 available on Dryad). Those haplotypes within var. ponderosa occupy environmental space that contains higher winter precipitation than haplotypes within var. scopulorum, in agreement with previous analyses (Norris et al. 2006; Shinneman et al. 2016; Fig. 6). Within var. ponderosa, H1 and H8 have the most overlap (D |$=$| 0.497), whereas H8 and H5 have very little overlap (D |$=$| 0.03). Within var. scopulorum, the most overlap is between H3 and H7 (D |$=$| 0.324) and H3 and H2 (D |$=$| 0.244) whereas H4 is the most isolated (D |$<$| 0.05 with all haplotypes within the variety). The niche space of H4, however, overlaps with the niche space of H1 (D |$=$| 0.126) and H8 (D |$=$| 0.260) of var. ponderosa (Fig. 6). Otherwise, haplotypes from var. ponderosa have very little niche overlap with haplotypes within var. scopulorum. H2, which has populations in New Mexico and southern California, has been placed within var. scopulorum (Potter et al. 2013) and the niche-overlap results support this (indeed, H2 overlaps with H3 [D |$=$| 0.244] more than other haplotypes). However, the low niche overlap between H2 and geographically proximate haplotypes in California associated with var. ponderosa, H1 (D |$=$| 0.023) and H5 (D |$=$| 0.014), could also be influenced by the lack of H2 presences predicted for southern California by the lineage distance modeling. Figure 6. Open in new tabDownload slide a) Niche of each variety and haplotype in climatic space defined by the two first axes of the PCA of the seven climate variables used in ENMs for western United States. The amount of variation explained by PC Axis 1 and 2 is 38.54% and 28.84%, respectively. Gray shading represents the density of occurrences of each variety and the lines indicate the outline of the niche space for each variety and haplotype. Inset is the contribution of the seven climate variables on the two axes of the PCA. b) The amount of niche overlap (D statistic) and |$P$|-values for niche similarity between each set of haplotypes (Supplementary Appendix S6 available on Dryad). Figure 6. Open in new tabDownload slide a) Niche of each variety and haplotype in climatic space defined by the two first axes of the PCA of the seven climate variables used in ENMs for western United States. The amount of variation explained by PC Axis 1 and 2 is 38.54% and 28.84%, respectively. Gray shading represents the density of occurrences of each variety and the lines indicate the outline of the niche space for each variety and haplotype. Inset is the contribution of the seven climate variables on the two axes of the PCA. b) The amount of niche overlap (D statistic) and |$P$|-values for niche similarity between each set of haplotypes (Supplementary Appendix S6 available on Dryad). Discussion Comparisons Among Projected ENM Distributions for Varieties Versus Haplotypes Phylogeographic structures defined by genetic variation may represent clearly defined and distinct niche spaces among lineages within a species that reflect the evolutionary capacity of a species to adapt to environmental change (Prasad and Potter 2017). Thus, modeling the ecological niches of intraspecific lineages may project distributions of each lineage that differ in important ways from projections for the species (Hällfors et al. 2016). For instance, using the fundamental niche of simulated conceptual species and subpopulations, Valladares et al. (2014) found more restricted distributions when modeling distinct subpopulations compared with modeling the entire species. In contrast, Pearman et al. (2010), Oney et al. (2013), and Valladares et al. (2014) found the opposite when modeling the realized niches of real species and populations within them. Valladares et al. (2014) attributed such differences to modeling the fundamental versus realized niches. However, our work suggests an alternative explanation is in order, as our realized niches of the haplotypes were generally more limiting than those of the varieties. Estimates derived from the collective representation of the haplotype probabilities for the recent period are also generally more geographically restricted and have higher average occurrence probabilities than those projected for the varieties (Figs. 2, 4, and 5), and these more restricted ranges are further amplified for most individual haplotypes under future climate change compared with the varieties (Supplementary Appendix S5 available on Dryad). Although this could be partly caused by the influence of haplotype presence locations that restrict estimates of the peripheral climate niche, we suggest that another plausible explanation is that the haplotype populations correspond to generally well defined and often minimally overlapping climate niches (Fig. 6) that are uniquely projected to decrease in the future. In contrast, models for the varieties capture a broader niche representing more potential combinations of climate conditions, a larger portion of which may persist into the future. It has been suggested that by not including intraspecific variation in ENMs, fitness-environment curves may be smoothed due to averaging across the broader niche, leading to more pessimistic forecasts of future distributions (Valladares et al. 2014). However, our variety models did not result in more unfavorable distributions (i.e., better delineated with higher probabilities) relative to the haplotype models. Regardless of intraspecific emphasis, by examining relative changes across models (based on probability estimates |$>$| 0.5), we were able to reveal potential trends that generally hold for both the variety and most corresponding haplotype models under climate change. For the Pacific variety (var. ponderosa), there is a general loss of area projected under climate change (by ~ 25%) and a modest increase in elevation (⁠|$<$| 200 m). Although individual haplotype responses were more variable, they generally supported trends of decreasing range size and increasing elevation projected by the variety model under future climates; however, the combined haplotype distributions did not (Figs. 2 and 3). Projected trends for the Rocky Mountain variety (var. scopulorum) were similar but more extreme (⁠|$>$| 50% decrease in area, and |$>$| 600 m elevation gain) than the Pacific variety, and were also generally supported by the more highly variable individual haplotype models, and to a lesser degree by the combined haplotype models (Figs. 3 and 5). These trends suggest that the individual haplotype models effectively captured inherent sub-specific variability in response, but that the combined haplotypes were less likely to decline overall because their ranges were already more restricted for the recent period (relative to ranges for the varieties), or because they were buoyed by a few favorable individual haplotype responses to climate change (e.g., projected distributions for H5, see Discussion below). We expected some predicted distribution loss in the future, especially given that extrapolation of ENMs fit with present-day climates and projected to novel climates tend to under-predict distributions (Maguire et al. 2016). This factor could be further enhanced when ENMs are fit with a narrower subset of occurrences, restricting potential future distributions. However, climate novelty is likely minimal in our models (⁠|$<$| 0.0003 following methods in Maguire et al. 2016), and the varying degrees of predicted area loss for haplotypes and varieties suggest potentially true signals in their projected responses to climate change. Indeed, the potentially unique responses of individual haplotypes to future climate within each ponderosa pine variety can be instructive, especially in terms of patterns of area loss and elevational shifts. Given likely long-term isolation and local evolution of each of the haplotypes within glacial periods, each may be expected to exhibit different responses to climate change (Potter et al. 2013; Shinneman et al. 2016). For example, within var. ponderosa, H5 does not lose a significant amount of area in the future under climate change (though probability values in estimated recent locations are projected to decline; Fig. 4), and it is the only haplotype projected to maintain a relatively stable mean elevation over time (Fig. 3). This may be due to ensemble climate projections that suggest relatively moderate warming and wet growing season conditions in the predominantly coastal and eastern Sierra Nevada range of this haplotype. Rehfeldt et al. (2014) also predicted relative geographic stability for var. ponderosa under future climate scenarios in this portion of its range. In contrast, within var. scopulorum, H6 is projected to lose most of its current range and shift significantly upwards (⁠|$>$| 1200 m) in mean elevation, most likely due to projections of hotter and drier growing seasons in the central and northeastern Rocky Mountains. The dramatic loss of H6 and var. scopulorum in the north (the lower elevations of the Black Hills and the hills and tablelands of south-eastern and central Montana) is in agreement with previous studies that modeled the niche distribution of the variety (Rehfeldt et al. 2014). The variety already occupies higher elevational areas in much of this part of its range, and there are large distances to the nearest analogous climate in the future, effectively restricting both elevational and latitudinal range shifts (Rehfeldt et al. 2014, Roberts and Hamann 2016). Haplotype 6 may also be vulnerable because it occupies a unique climate niche space (e.g., occupies higher values of growing season precipitation, precipitation ratio, and growing degree days below 0|$^{\circ}$|C than other haplotypes) that is on the edge of the total niche space for var. scopulorum (Fig. 6 and Supplementary Appendix S6 available on Dryad). Relevance to Adaptation and Phenotypic Plasticity To successfully adapt to climate change, plant species must either migrate to track their ecological niches or increase their tolerance to new climate conditions in situ. This latter strategy requires either expression of phenotypic plasticity, or selection on existing phenotypes defined by existing genetic material or on novel phenotypes resulting from favorable mutations (Williams et al. 2008; Chevin et al. 2010). Our haplotype ENMs operated under an assumption of adaptation to specific, realized climate niches; they did not directly incorporate plasticity or genetic adaptation under changing climate, and thus did not capture the ability of each haplotype to expand into or occupy other parts of the species’ fundamental niche. Moreover, we caution that the haplotypes, which are classified based on identification of neutral genetic markers, are not necessarily effective surrogates for uniquely adapted populations (e.g., Holderegger et al. 2006). However, although these mtDNA defined haplotypes do not necessarily have adaptive significance, they correspond to populations that have been influenced by long-term biogeographical processes. Previous studies of ponderosa pine have shown the existence of adaptive phenotypic variation within populations (e.g. Conkle 1973; Rehfeldt 1986a, 1986b, 1990; Namkoong and Conkle 1976; Sorensen and Weber 1994; Sorensen et al. 2001; Keller et al. 2004), and populations with greater phenotypic plasticity may be more tolerant of changing climates (Chevin et al. 2010). There are areas within the distribution of ponderosa pine with high genetic variation (based on nuclear DNA, see Potter et al. 2015), and greater genetic variation can allow populations to adapt to climate change via natural selection (Sork et al. 2010). In addition, high gene flow has been detected within each ponderosa pine variety, occurring across complex structures of genetic variation (Potter et al. 2015). Gene flow between even the most genetically structured populations can reduce population differentiation and relax migrational constraints (Potter et al. 2015). Given the potential for plasticity or genetically adaptive responses of ponderosa pine to changing climate, as well as shifts in biotic interactions, population-level responses may be best captured by developing ENMs for the varieties. Indeed, when reduced to two dimensions via the niche-overlap analysis, the collective ponderosa pine haplotype climate niche space does not cover the entire niche space of each variety (Fig. 6), potentially suggesting that the modeled niches of the varieties could indirectly represent more plasticity or adaptive potential for the species than the individual haplotypes combined. The ENMs of varieties may also be indirectly capturing intraspecific variation, and therefore, the plasticity of the variety, reflected in the projected broader geographic distributions and higher probabilities under climate change compared with most haplotypes. This highlights the complementary value of modeling at both the variety and haplotype levels, in relation to phenotypic plasticity and adaptation respectively, and providing a better understanding of how variability in intraspecific response might affect the future niche distribution of the species as a whole (Pauls et al. 2013; Marcer et al. 2016). Dispersal, Migration, and Competition Even though optimal climate for ponderosa pine is generally predicted to shift to higher elevations in the future, existing populations may not fill these areas because of dispersal limitations and migration lags (Davis 1989). Thus, rates of migration on leading edges may be slower than rates of habitat loss on trailing edges of existing populations (Campbell and Shinneman 2017). Yet, although ponderosa pine has limited optimal seed dispersal distance for regeneration (15–30 m; Oliver and Ryker 1990; Boyden et al. 2005; Latta et al. 1998), it is capable of low frequency, long-distance dispersal of 75–200 km (Lesser and Jackson 2012, 2013). Indeed, this likely explains its rapid expansion and northward migration from refugial areas during the Last Glacial Maximum, as supported by fossil evidence (Norris et al. 2016) and by a pattern of decreasing unique and rare alleles with increasing latitude found in contemporary populations (Potter et al. 2015). These paleogeographic range shifts may suggest potential for a similar response to future climate change. However, requirements of future migrations may be different from the past due to novel climate conditions, different refugia locations, and new migrational constraints imposed by human activity (Norris et al. 2016; Roberts and Hamann 2016). In addition, because realized climate niches (from which the ENMs were built) may be truncated due to biotic interactions; shifting interspecific competition under climate change could also affect future distributions (Wisz et al. 2013). Potential Genetic Variation and Diversity Under Future Climate Genetic variation and diversity are likely to be affected differently for ponderosa pine populations under future and potentially novel climate conditions compared with paleoclimates. For instance, based on fossil and genetic evidence (Potter et al. 2013), the contact zone between the two varieties in central Montana is relatively recent, established within the previous 1500 years, and likely controlled by a geographic shift from summer- to winter-dominant precipitation (Norris et al. 2006). Chloroplast and nuclear DNA patterns suggest introgression from west to east in this region (Latta and Mitton 1999; Potter et al. 2015), but there was no exchange of mtDNA (Latta and Mitton 1999; Johansen and Latta 2003; Potter et al. 2013). Potential genetic exchange between haplotypes and varieties in this location could be further restricted under future climate, based on the lack of overlap projected in our models for either the varieties or haplotypes. In addition, some haplotypes could become increasingly isolated under warmer and drier climates (e.g., H4 and H7). Thus, some genetic diversity may be lost, especially if isolated lineages have less suitable habitat and experience even greater geographic isolation under future climates. This underscores the need to identify potential future refugia across intraspecific levels that may offer the best chance for maintenance of genetic diversity and survival of the species under climate change (e.g., Keppel et al. 2012). Conclusions Our results highlight the potential value of modeling intraspecific response to climate change and, more specifically, are informative for future management and conservation of ponderosa pine, an economically, and ecologically important tree species. Ponderosa pine haplotype groups may represent evolutionary lineages adapted to different environmental conditions. Thus, by identifying lineages most vulnerable to climate change and detecting potential climatically suitable areas in the future, our results suggest at least three important considerations for successful conservation. First, it is significant that each haplotype occupies a statistically different environmental space within the range of ponderosa pine, as this could be considered evidence for differential adaptation to diverse environmental conditions and, thus, that the haplotypes may be evolutionary units worthy of distinct conservation attention. However, there is need for further study on adaptive genetic variation in ponderosa pine that could help to identify evolutionarily significant units for conservation, including defining genetic groups differently (i.e., based on nuclear DNA variation; Potter et al. 2015). Second, the results indicate that some haplotype groups, specifically H8 in var. ponderosa and H6 in var. scopulorum, could lose well more than half of their predicted distribution before 2090. These haplotypes may require special efforts to conserve their existing genetic variation, using both ex situ and in situ strategies (Potter et al. 2017). Peripheral populations of isolated lineages, such as H4 and H7, may also be at high risk under climate change, and due to their lower genetic diversity but greater genetic differentiation compared with populations in the core range of each ponderosa pine variety (Potter et al. 2015), may be ideal candidates for conservation action. Finally, the results of this study could help guide assisted species and population migration efforts (Dumroese et al. 2015; Aitken and Bemmels 2016), by identifying locations in which vulnerable ponderosa pine lineages may be best adapted in the future. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.jj264. Funding This work was supported in part by Cost Shares Agreement [14-CS-11330110-042 and 15-CS-11330110-067 to K.M.P.] between the Southern Research Station of the USDA Forest Service and North Carolina State University. Funding for developing the climate niche models and niche-overlap analysis presented in this article was provided by an Interagency Agreement [L12PG00219 to D.J.S.] between the U.S. Bureau of Land Management and the U.S. Geological Survey. Acknowledgments Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government. We would like to thank Dan Rosauer, Nicholas Matzke, and Brittany Barker for their very helpful reviews of this article. We would also like to acknowledge Robert Means, whose efforts helped to initiate this research, and who passed away before the work was completed. References Aitken S.N, Bemmels J.B. 2016 . Time to get moving: assisted gene flow of forest trees. Evol. Appl. 9 : 271 – 290 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Bálint M., Domisch S., Engelhardt C.H.M., Haase P., Lehrian S., Sauer J., Theissinger K., Pauls S.U., Nowak C. 2011 . Cryptic biodiversity loss linked to global climate change. Nat. Clim. Change 1 : 313 – 318 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Benito Garzón M., Alía R., Robson T.M., Zavala M.A. 2011 . Intra-specific variability and plasticity influence potential tree species distributions under climate change: intra-specific variability and plasticity. Global Ecol. Biogeogr. 20 : 766 – 778 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Boyden S., Binkley D., Shepperd W. 2005 . Spatial and temporal patterns in structure, regeneration, and mortality of an old-growth ponderosa pine forest in the Colorado Front Range. Forest Ecol. Manag. 219 : 43 – 55 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Brier G.W. 1950 . Verification of forecasts expressed in terms of probability. Mon. Weather Rev. 78 : 1 – 3 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Broennimann O., Fitzpatrick M.C., Pearman P.B., Petitpierre B., Pellissier L., Yoccoz N.G., Thuiller W., Fortin, M.-J., Randin C., Zimmermann N.E., Graham C.H., Guisan A. 2012 . Measuring ecological niche overlap from occurrence and spatial environmental data: Measuring niche overlap. Global Ecol. Biogeogr. 21 : 481 – 497 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Broennimann O., Di Cola V., Guisan A. 2016 . ecospat: spatial ecology miscellaneous methods. R package version 2.1.1. Available from: URL https://CRAN.R-project.org/package=ecospat. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Campbell J.L., Shinneman D.J. 2017 . Potential influence of wildfire in modulating climate-induced forest redistribution in a central Rocky Mountain landscape. Ecol. Process. 6 : 7 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Chevin L.-M., Lande R., Mace G.M. 2010 . Adaptation, plasticity, and extinction in a changing environment: towards a predictive theory. PLoS Biol. 8 : e1000357 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Conkle M.T. 1973 . Growth data for 29 years from the California elevational transect study of ponderosa pine. Forest Sci. 19 : 31 – 39 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crookston N.L., Rehfeldt G.E. 2008 . Climate estimates and plant-climate relationships. Moscow (Idaho) : USDA Forest Service, Rocky Mountain Research Station. Available from: URL http://forest.moscowfsl.wsu.edu/climate/. Google Scholar D’Amen M., Zimmermann N.E., Pearman P.B. 2013 . Conservation of phylogeographic lineages under climate change: African mammals and global warming. Global Ecol. Biogeogr. 22 : 93 – 104 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Davis M.B. 1989 . Lags in vegetation response to greenhouse warming. Climatic Change 15 : 75 – 82 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Dumroese R.K., Williams M.I., Stanturf J.A., St. Clair J.B. 2015 . Considerations for restoring temperate forests of tomorrow: forest restoration, assisted migration, and bioengineering. New Forests 46 : 947 – 964 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Espíndola A., Pellissier L., Maiorano L., Hordijk W., Guisan A., Alvarez N. 2012 . Predicting present and future intra-specific genetic structure through niche hindcasting across 24 millennia: hindcasting-based phylogeography. Ecol. Lett. 15 : 649 – 657 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hällfors M.H., Liao J., Dzurisin J., Grundel R., Hyvärinen M., Towle K., Wu G.C., Hellmann J.J. 2016 . Addressing potential local adaptation in species distribution models: implications for conservation under climate change. Ecol. Appl. 26 : 1154 – 1169 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Holderegger R., Kamm U., Gugerli F. 2006 . Adaptive vs. neutral genetic diversity: implications for landscape genetics. Landsc. Ecol. 21 : 797 – 807 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Jaramillo-Correa J.P., Beaulieu J., Khasa D.P., Bousquet J. 2009 . Inferring the past from the present phylogeographic structure of North American forest trees: seeing the forest for the genes. Can. J. Forest Res. 39 : 286 – 307 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Johansen A.D., Latta R.G. 2003 . Mitochondrial haplotype distribution, seed dispersal and patterns of postglacial expansion of ponderosa pine. Mol. Ecol. 12 : 293 – 298 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Keller E.A., Anekonda T.S., Smith B.N., Hansen L.D., Clair J.B.S., Criddle R.S. 2004 . Stress and respiration traits differ among four geographically distinct Pinus ponderosa seed sources. Thermochim. Acta 422 : 69 – 74 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Keppel G., Van Niel K.P., Wardell-Johnson G.W., Yates C.J., Byrne M. Mucina L., Schut A.G.T., Hopper S.D., Franklin S.E. 2012 . Refugia: identifying and understanding safe havens for biodiversity under climate change. Global Ecol. Biogeogr. 21 : 393 – 404 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lascoux M., Palmé A.E., Cheddadi R., Latta R.G. 2004 . Impact of Ice Ages on the genetic structure of trees and shrubs. Philos. Trans. R. Soc. Lond. B Biol. Sci. 359 : 197 – 207 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Latta R.G., Linhart Y.B., Fleck D., Elliot M. 1998 . Direct and indirect estimates of seed versus pollen movement within a population of ponderosa pine. Evolution 52 : 61 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Latta R.G., Mitton J.B. 1999 . Historical separation and present gene flow through a zone of secondary contact in ponderosa pine. Evolution 53 : 769 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lesser M.R., Jackson S.T. 2012 . Making a stand: five centuries of population growth in colonizing populations of Pinus ponderosa. Ecology 93 : 1071 – 1081 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lesser M.R., Jackson S.T. 2013 . Contributions of long-distance dispersal to population growth in colonising Pinus ponderosa populations. Ecol. Lett. 16 : 380 – 389 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Liaw A., Wiener M. 2002 . Classification and regression by randomForest. R News 2 : 18 – 22 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Little E.L. 1971 . Atlas of United States Trees, vol. 1 . Conifers and important hardwoods. Washington (DC) : USDA Forest Service. Miscellaneous Publication 1146. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Maguire K.C., Nieto-Lugilde D., Blois J.L., Fitzpatrick M.C., Williams J.W., Ferrier S., Lorenz D.J. 2016 . Controlled comparison of species- and community-level models across novel climates and communities. Proc. R. Soc. Lond. B Biol. Sci. 283 : 20152817 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Marcer A., Méndez-Vigo B., Alonso-Blanco C., Picó F.X. 2016 . Tackling intraspecific genetic structure in distribution models better reflects species geographical range. Ecol. Evol. 6 : 2084 – 2097 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC McCune B., and Mefford J. 2009 . HyperNiche. Nonparametric multiplicative habitat modeling, Version 2.28. Gleneden Beach (OR) : MjM Software. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Milborrow S. 2016 . Earth: multivariate adaptive regression splines. R package version 4.4.4. Available from: URL https://CRAN.R-project.org/package=earth. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Namkoong G., Conkle M.T. 1976 . Time trends in genetic control of height growth in ponderosa pine. Forest Sci. 22 : 2 – 12 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Neale D.B., Sederoff R.R. 1989 . Paternal inheritance of chloroplast DNA and maternal inheritance of mitochondrial DAN in loblolly pine. Theor. Appl. Genet. 77 : 212 – 216 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Niebling C.R., Conkle M.T. 1990 . Diversity of Washoe pine and comparisons with allozymes of ponderosa pine races. Can. J. Forest Res. 20 : 298 – 308 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Norris J.R., Jackson S.T., Betancourt J.L. 2006 . Classification tree and minimum-volume ellipsoid analyses of the distribution of ponderosa pine in the western USA. J. Biogeogr. 33 : 342 – 360 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Norris J.R., Betancourt J.L., Jackson S.T. 2016 . Late Holocene expansion of ponderosa pine (Pinus ponderosa) in the Central Rocky Mountains, USA. J. Biogeogr. 43 : 778 – 790 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Oliver W.W., Ryker R.A. 1990 . Ponderosa pine. In: Burns R.M., Honkala B.H., editors. Silvics of North America: 1. Conifers, vol 1 . Agricultural Handbook 654. Washington : U.S. Department of Agriculture Forest Service. p. 413 – 424 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Oney B., Reineking B., O’Neill G., Kreyling J. 2013 . Intraspecific variation buffers projected climate change impacts on Pinus contorta. Ecol. Evol. 3 : 437 – 449 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pauls S.U., Nowak C., Bálint M., Pfenninger M. 2013 . The impact of global climate change on genetic diversity within populations and species. Mol. Ecol. 22 : 925 – 946 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pearman P.B., D’Amen M., Graham C.H., Thuiller W., Zimmermann N.E. 2010 . Within-taxon niche structure: niche conservatism, divergence and predicted effects of climate change. Ecography 33 : 990 – 1003 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Petit R.J., Kremer A., Wagner B. 1993 . Finite island model for organelle and nuclear genes in plants. Heredity 71 : 630 – 641 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pfenninger M., Nowak C., Magnin F. 2007 . Intraspecific range dynamics and niche evolution in Candidula land snail species. Biol. J. Linn. Soc. 90 : 303 – 317 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Potter K.M., Hipkins V.D., Mahalovich M.F., Means R.E. 2013 . Mitochondrial DNA haplotype distribution patterns in Pinus ponderosa (Pinaceae): range-wide evolutionary history and implications for conservation. Am. J. Bot. 100 : 1562 – 1579 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Potter K.M., Hipkins V.D., Mahalovich M.F., Means R.E. 2015 . Nuclear genetic variation across the range of ponderosa pine (Pinus ponderosa): phylogeographic, taxonomic and conservation implications. Tree Genet. Genomes 11 : 38 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Potter K.M., Jetton R.M., Bower A., Jacobs D.F., Man G., Hipkins V.D., Westwood M. 2017 . Banking on the future: progress, challenges and opportunities for the genetic conservation of forest trees. New Forests 48 : 153 – 180 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Prasad A.M., Potter K.M. 2017 . Macro-scale assessment of demographic and environmental variation within genetically derived evolutionary lineages of eastern hemlock (Tsuga canadensis), an imperiled conifer of the eastern United States. Biodivers. Conserv. 10.1007/s10531-017-1354-4 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref R Core Team 2014 . R: a language and environment for statistical computing. Vienna (Austria) : R Foundation for Statistical Computing. Available from: URL http://www.R-project.org/. Rehfeldt G.E. 1986b . Adaptive variation in Pinus ponderosa from Intermountain regions. II. Middle Columbia River System. Research Paper INT-RP-373. Ogden, UT : U.S. Department of Agriculture, Forest Service, Intermountain Research Station. p9 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E. 1986b . Adaptive variation in Pinus ponderosa from intermountain regions. II. Middle Columbia river system. Research Paper INT-RP-373. Ogden, UT : U.S. Department of Agriculture, Forest Service, Intermountain Research Station. p9 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E. 1990 . Genetic differentiation among populations of Pinus ponderosa from the Upper Colorado River Basin. Bot. Gaz. 151 : 125 – 137 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E. 2006 . A spline model of climate for the Western United States. General Technical Report RMRS-GTR-165. Fort Collins, CO : U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station. p.21 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E., Jaquish B.C., López-Upton J., Sáenz-Romero C., St Clair J.B., Leites L.P., Joyce D.G. 2014 . Comparative genetic responses to climate for the varieties of Pinus ponderosa and Pseudotsuga menziesii: realized climate niches. Forest Ecol. Manag. 324 : 126 – 137 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rosauer D.F., Catullo R.A., VanDerWal J., Moussalli A., Moritz C. 2015 . Lineage range estimation method reveals fine-scale endemism linked to Pleistocene stability in Australian rainforest herpetofauna. PLoS One 10 : e0126274 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Roberts D.R., Hamann A. 2016 . Climate refugia and migration requirements in complex landscapes. Ecography 39 : 1238 – 1246 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rudloff E., von Lapp M.S. 1992 . Chemosystematic studies in the genus Pinus. VII. The leaf oil terpene composition of ponderosa pine, Pinus ponderosa. Can. J. Bot. 70 : 374 – 378 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Shinneman D.J., Means R.E., Potter K.M., Hipkins V.D. 2016 . Exploring climate niches of ponderosa pine (Pinus ponderosa Douglas ex Lawson) haplotypes in the western United States: implications for evolutionary history and conservation. PLoS One 11 : e0151811 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Schoener T.W. 1970 . Nonsynchronous spatial overlap of lizards in patchy habitats. Ecology 51 : 408 – 418 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sorensen F.C., Mandel N.L., Aagaard J.E. 2001 . Role of selection versus historical isolation in racial differentiation of ponderosa pine in southern Oregon: an investigation of alternative hypotheses. Can. J. Forest Res. 31 : 1127 – 1139 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sorensen F.C., Weber J.C. 1994 . Genetic variation and seed transfer guidelines for ponderosa pine in the Ochoco and Malheur National Forests of central Oregon. USDA Forest Service Research Paper PNW-RP-472, Washington DC , USA . 24 p. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sork V.L., Davis F.W., Westfall R., Flint A., Ikegami M., Wang H., Grivet D. 2010 . Gene movement and genetic association with regional climate gradients in California valley oak (Quercus lobata Née) in the face of climate change. Mol. Ecol. 19 : 3806 – 3823 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Valladares F., Matesanz S., Guilhaumon F., Araújo M.B., Balaguer L., Benito-Garzón M., Cornwell W., Gianoli E., van Kleunen M., Naya D.E., Nicotra A.B., Poorter H., Zavala M.A. 2014 . The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change. Ecol. Lett. 17 : 1351 – 1364 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wang C.-W. 1977 . Genetics of ponderosa pine. Washington (DC) : U.S. Department of Agriculture, Forest Service . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wells O.O. 1964 . Geographic variation in Ponderosa Pine. Silvae Genetica 13 : 89 – 103 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Williams S.E., Shoo L.P., Isaac J.L., Hoffmann A.A., Langham G. 2008 . Towards an integrated framework for assessing the vulnerability of species to climate change. PLoS Biol. 6 : e325 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wisz M.S., Pottier J., Kissling W.D., Pellissier L., Lenoir J., Damgaard C.F., Dormann C.F., Forchhammer M.C., Grytnes J.-A., Guisan A., Heikkinen R.K., Høye T.T., Luoto M., Maiorano L., Nilssson M.-C., Normand S., Öckinger E., Schmidt N.M., Termansen M., Timmerman A., Wardle D.A., Aastrup P., Svenning J.-C. 2013 . The role of biotic interactions in shaping distributions and realized assemblages of species: implications for species distribution modelling. Biol. Rev. 88 : 15 – 30 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wofford A.M., Finch K., Bigott A., Willyard A. 2014 . A set of plastid loci for use in multiplex fragment length genotyping for intraspecific variation in Pinus (Pinaceae). Appl. Plant Sci. 2 : 1400002 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Yannic G., Pellissier L., Ortego J., Lecomte N., Couturier S., Cuyler C., Dussault C., Hundertmark K.J., Irvine R.J., Jenkins D.A., Kolpashikov L., Mager K., Musiani M., Parker K.L., Røed K.H., Sipko T., Þórisson S.G., Weckworth B.V., Guisan A., Bernatchez L., Côté S.D. 2014 . Genetic diversity in caribou linked to past and future climate change. Nat. Clim. Change 4 : 132 – 137 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Published by Oxford University Press, on behalf of the Society of Systematic Biologists 2018. This work is written by (a) US Government employee(s) and is in the public domain in the US. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) Published by Oxford University Press, on behalf of the Society of Systematic Biologists 2018. This work is written by (a) US Government employee(s) and is in the public domain in the US. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Systematic Biology Oxford University Press http://www.deepdyve.com/lp/oxford-university-press/intraspecific-niche-models-for-ponderosa-pine-pinus-ponderosa-suggest-eLi6MZavM7

Loading next page...

References (121)

F. Valladares, S. Matesanz, F. Guilhaumon, M.B. Araújo, L. Balaguer, M. Benito-Garzón, W. Cornwell, E. Gianoli, M. van Kleunen, D.E. Naya, A.B. Nicotra, H. Poorter, M.A. Zavala (2014)
Ecol. Lett.
, 17
Espíndola (2012)
Predicting present and future intra-specific genetic structure through niche hindcasting across 24 millennia: hindcasting-based phylogeography.
K. Potter, V. Hipkins, M. Mahalovich, R. Means (2013)
Mitochondrial DNA haplotype distribution patterns in Pinus ponderosa (Pinaceae): range-wide evolutionary history and implications for conservation.
American journal of botany, 100 8
S.N Aitken, J.B. Bemmels (2016)
Evol. Appl.
, 9
D.R. Roberts, A. Hamann (2016)
Ecography
, 39
M. D’Amen, N. Zimmermann, P. Pearman (2013)
Conservation of phylogeographic lineages under climate change
Global Ecology and Biogeography, 22
K. Potter, R. Jetton, A. Bower, D. Jacobs, G. Man, V. Hipkins, Murphy Westwood (2017)
Banking on the future: progress, challenges and opportunities for the genetic conservation of forest trees
New Forests, 48
J. Jaramillo‐Correa, J. Beaulieu, D. Khasa, J. Bousquet (2009)
Inferring the past from the present phylogeographic structure of North American forest trees : seeing the forest for the genes
Canadian Journal of Forest Research, 39
G. Rehfeldt (2006)
A spline model of climate for the Western United States
, 165
A. Prasad, K. Potter (2017)
Macro-scale assessment of demographic and environmental variation within genetically derived evolutionary lineages of eastern hemlock (Tsuga canadensis), an imperiled conifer of the eastern United States
Biodiversity and Conservation, 26
R. Holderegger, U. Kamm, F. Gugerli (2006)
Adaptive vs. neutral genetic diversity: implications for landscape genetics
Landscape Ecology, 21
A.M. Wofford, K. Finch, A. Bigott, A. Willyard (2014)
Appl. Plant Sci.
, 2
R. Holderegger, U. Kamm, F. Gugerli (2006)
Landsc. Ecol.
, 21
M.S. Wisz, J. Pottier, W.D. Kissling, L. Pellissier, J. Lenoir, C.F. Damgaard, C.F. Dormann, M.C. Forchhammer, J.-A. Grytnes, A. Guisan, R.K. Heikkinen, T.T. Høye, M. Luoto, L. Maiorano, M.-C. Nilssson, S. Normand, E. Öckinger, N.M. Schmidt, M. Termansen, A. Timmerman, D.A. Wardle, P. Aastrup, J.-C. Svenning (2013)
Biol. Rev.
, 88
V.L. Sork, F.W. Davis, R. Westfall, A. Flint, M. Ikegami, H. Wang, D. Grivet (2010)
Mol. Ecol.
, 19
M. Wisz, J. Pottier, Werner Kissling, L. Pellissier, J. Lenoir, C. Damgaard, C. Dormann, M. Forchhammer, J. Grytnes, A. Guisan, R. Heikkinen, T. Høye, I. Kühn, M. Luoto, L. Maiorano, M. Nilsson, S. Normand, Erik Öckinger, N. Schmidt, M. Termansen, A. Timmermann, D. Wardle, P. Aastrup, J. Svenning (2012)
The role of biotic interactions in shaping distributions and realised assemblages of species: implications for species distribution modelling
Biological Reviews of the Cambridge Philosophical Society, 88
(1950)
Mon. Weather Rev.
, 78
(2011)
Intraspecific variability and plasticity influence potential tree species distributions under climate change: intra-specific variability and plasticity
J.L. Campbell, D.J. Shinneman (2017)
Ecol. Process.
, 6
M. Conkle (1973)
Growth Data for 29 Years from the California Elevational Transect Study of Ponderosa Pine
Forest Science, 19
(1990)
River Basin. Bot. Gaz.
, 151
B. McCune, J. Mefford (2009)
Nonparametric multiplicative habitat modeling, Version 2.28.
K.C. Maguire, D. Nieto-Lugilde, J.L. Blois, M.C. Fitzpatrick, J.W. Williams, S. Ferrier, D.J. Lorenz (2016)
Proc. R. Soc. Lond. B Biol. Sci.
, 283
S. Aitken, Jordan Bemmels (2015)
Time to get moving: assisted gene flow of forest trees
Evolutionary Applications, 9
(2006)
General Technical Report RMRS-GTR-165.
J.R. Norris, J.L. Betancourt, S.T. Jackson (2016)
J. Biogeogr.
, 43
O. Broennimann, V. Di Cola, A. Guisan (2016)
R package version 2.1.1.
G. Keppel, K. Niel, G. Wardell-Johnson, C. Yates, M. Byrne, L. Mucina, A. Schut, S. Hopper, S. Franklin (2012)
Refugia: identifying and understanding safe havens for biodiversity under climate change
Global Ecology and Biogeography, 21
S. Pauls, C. Nowak, Miklós Bálint, M. Pfenninger (2013)
The impact of global climate change on genetic diversity within populations and species
Molecular Ecology, 22
G.E. Rehfeldt, B.C. Jaquish, J. López-Upton, C. Sáenz-Romero, Clair J.B. St, L.P. Leites, D.G. Joyce (2014)
Forest Ecol. Manag.
, 324
A. Espíndola, L. Pellissier, L. Maiorano, W. Hordijk, A. Guisan, N. Alvarez (2012)
Predicting present and future intra-specific genetic structure through niche hindcasting across 24 millennia.
Ecology letters, 15 7
M. Lascoux, A.E. Palmé, R. Cheddadi, R.G. Latta (2004)
Philos. Trans. R. Soc. Lond. B Biol. Sci.
, 359
M. Lesser, S. Jackson (2012)
Making a stand: five centuries of population growth in colonizing populations of Pinus ponderosa.
Ecology, 93 5
A. Johansen, R. Latta (2002)
Mitochondrial haplotype distribution, seed dispersal and patterns of postglacial expansion of ponderosa pine
Molecular Ecology, 12
G. Rehfeldt, B. Jaquish, J. López-Upton, C. Sáenz-Romero, J. Clair, L. Leites, D. Joyce (2014)
Comparative genetic responses to climate for the varieties of Pinus ponderosa and Pseudotsuga menziesii: Realized climate niches
Forest Ecology and Management, 324
G. Namkoong, M.T. Conkle (1976)
Forest Sci.
, 22
T. Schoener (1970)
Nonsynchronous Spatial Overlap of Lizards in Patchy Habitats
Ecology, 51
Milborrow (2016)
Earth: multivariate adaptive regression splines.
M. Hällfors, Jishan Liao, Jason Dzurisin, R. Grundel, M. Hyvärinen, Kevin Towle, Grace Wu, J. Hellmann (2016)
Addressing potential local adaptation in species distribution models: implications for conservation under climate change.
Ecological applications : a publication of the Ecological Society of America, 26 4
(1971)
Conifers and important hardwoods.
, vol. 1
(1964)
Silvae Genetica
, 13
(1989)
Climatic Change
, 15
Austin Wofford, Kristen Finch, A. Bigott, A. Willyard (2014)
A set of plastid loci for use in multiplex fragment length genotyping for intraspecific variation in Pinus (Pinaceae)1
Applications in Plant Sciences, 2
(1977)
Washington (DC)
M.H. Hällfors, J. Liao, J. Dzurisin, R. Grundel, M. Hyvärinen, K. Towle, G.C. Wu, J.J. Hellmann (2016)
Ecol. Appl.
, 26
E.A. Keller, T.S. Anekonda, B.N. Smith, L.D. Hansen, J.B.S. Clair, R.S. Criddle (2004)
Acta
, 422
K.M. Potter, V.D. Hipkins, M.F. Mahalovich, R.E. Means (2015)
Genomes
, 11
J. Norris, J. Betancourt, S. Jackson (2016)
Late Holocene expansion of ponderosa pine (Pinus ponderosa) in the Central Rocky Mountains, USA
Journal of Biogeography, 43
A. Marcer, Bel Endez-Vigo, C. Alonso‐Blanco, F Pic, Correspondence Pic (2016)
Tackling intraspecific genetic structure in distribution models better reflects species geographical range
Ecology and Evolution, 6
C. Niebling, M. Conkle (1990)
Diversity of Washoe pine and comparisons with allozymes of ponderosa pine races
Canadian Journal of Forest Research, 20
Chi-Wu Wang (1977)
Genetics of ponderosa pine
G. Brier (1950)
VERIFICATION OF FORECASTS EXPRESSED IN TERMS OF PROBABILITY
Monthly Weather Review, 78
(1970)
Ecology
, 51
R.G. Latta, J.B. Mitton (1999)
Evolution
, 53
R. Team (2014)
R: A language and environment for statistical computing.
MSOR connections, 1
G. Yannic, L. Pellissier, J. Ortego, N. Lecomte, S. Couturier, C. Cuyler, C. Dussault, K.J. Hundertmark, R.J. Irvine, D.A. Jenkins, L. Kolpashikov, K. Mager, M. Musiani, K.L. Parker, K.H. Røed, T. Sipko, S.G. Þórisson, B.V. Weckworth, A. Guisan, L. Bernatchez, S.D. Côté (2014)
Nat. Clim. Change
, 4
F.C. Sorensen, N.L. Mandel, J.E. Aagaard (2001)
Can. J. Forest Res.
, 31
J. Campbell, D. Shinneman (2017)
Potential influence of wildfire in modulating climate-induced forest redistribution in a central Rocky Mountain landscape
Ecological Processes, 6
D’Amen (2013)
Conservation of phylogeographic lineages under climate change: African mammals and global warming.
R. Petit, A. Kremer, D. Wagner (1993)
Finite island model for organelle and nuclear genes in plants
Heredity, 71
M. Pfenninger, C. Nowak, F. Magnin (2007)
Biol. J. Linn. Soc.
, 90
Broennimann (2012)
Measuring ecological niche overlap from occurrence and spatial environmental data: Measuring niche overlap.
D. Shinneman, R. Means, K. Potter, V. Hipkins (2016)
Exploring Climate Niches of Ponderosa Pine (Pinus ponderosa Douglas ex Lawson) Haplotypes in the Western United States: Implications for Evolutionary History and Conservation
PLoS ONE, 11
F.C. Sorensen, J.C. Weber (1994)
USDA Forest Service Research Paper PNW-RP-472,
McCune (2009)
HyperNiche.
G. Yannic (2014)
Genetic diversity in caribou linked to past and future climate change
Nature Climate Change, 4
G. Keppel, K.P. Van Niel, G.W. Wardell-Johnson, C.J. Yates, M. Byrne, L. Mucina, A.G.T. Schut, S.D. Hopper, S.E. Franklin (2012)
Global Ecol. Biogeogr.
, 21
Rehfeldt (1986)
Adaptive variation in Pinus ponderosa from intermountain regions. II.
M. Pfenninger, C. Nowak, F. Magnin (2007)
Intraspecific range dynamics and niche evolution in Candidula land snail species
Biological Journal of The Linnean Society, 90
Emily Keller, T. Anekonda, Bruce Smith, L. Hansen, J. Clair, R. Criddle (2004)
Stress and respiration traits differ among four geographically distinct Pinus ponderosa seed sources
Thermochimica Acta, 422
G. Rehfeldt (1990)
Genetic Differentiation among Populations of Pinus ponderosa from the Upper Colorado River Basin
Botanical Gazette, 151
(2016)
R package version 4.4.4.
K.M. Potter, R.M. Jetton, A. Bower, D.F. Jacobs, G. Man, V.D. Hipkins, M. Westwood (2017)
New Forests
, 48
F. Sorensen, J. Weber (1994)
Genetic variation and seed transfer guidelines for ponderosa pine in the Ochoco and Malheur National Forests of central Oregon. Forest Service research paper
S. Boyden, D. Binkley, W. Shepperd (2005)
Spatial and temporal patterns in structure, regeneration, and mortality of an old-growth ponderosa pine forest in the Colorado Front Range
Forest Ecology and Management, 219
B. Oney, B. Reineking, G. O’Neill, J. Kreyling (2013)
Ecol. Evol.
, 3
G. Rehfeldt (1986)
Adaptive variation in Pinus ponderosa from Intermountain Regions. I. Snake and Salmon River Basins
Forest Science, 32
K.M. Potter, V.D. Hipkins, M.F. Mahalovich, R.E. Means (2013)
Am. J. Bot.
, 100
A.M. Prasad, K.M. Potter (2017)
Biodivers. Conserv.
A. Liaw, M. Wiener (2002)
R News
, 2
Luis-Miguel Chevin, R. Lande, G. Mace (2010)
Adaptation, Plasticity, and Extinction in a Changing Environment: Towards a Predictive Theory
PLoS Biology, 8
Kaitlin Maguire, D. Nieto-Lugilde, J. Blois, M. Fitzpatrick, John Williams, S. Ferrier, D. Lorenz (2016)
Controlled comparison of species- and community-level models across novel climates and communities
Proceedings of the Royal Society B: Biological Sciences, 283
K. Potter, V. Hipkins, M. Mahalovich, R. Means (2015)
Nuclear genetic variation across the range of ponderosa pine (Pinus ponderosa): Phylogeographic, taxonomic and conservation implications
Tree Genetics & Genomes, 11
M. Lesser, S. Jackson (2013)
Contributions of long-distance dispersal to population growth in colonising Pinus ponderosa populations.
Ecology letters, 16 3
S. Williams, L. Shoo, J. Isaac, A. Hoffmann, G. Langham (2008)
Towards an Integrated Framework for Assessing the Vulnerability of Species to Climate Change
PLoS Biology, 6
E. Little (1971)
Atlas of United States trees.
G. Namkoong, M. Conkle (1976)
Time trends in genetic control of height growth in Ponderosa Pine.
Forest Science, 22
E. Rudloff, M.S. von Lapp (1992)
Can. J. Bot.
, 70
R.J. Petit, A. Kremer, B. Wagner (1993)
Heredity
, 71
F. Sorensen, J. Weber (1994)
Genetic variation and seed transfer guidelines for ponderosa pine in the Ochoco and Malheur National Forests of central Oregon.
, 468
D. Roberts, D. Roberts, A. Hamann (2016)
Climate refugia and migration requirements in complex landscapes
Ecography, 39
Miklós Bálint, S. Domisch, C. Engelhardt, P. Haase, Stephanie Lehrian, Jan Sauer, K. Theissinger, S. Pauls, C. Nowak (2011)
Cryptic biodiversity loss linked to global climate change
Nature Climate Change, 1
(1986b)
Middle Columbia river system.
V. Sork, F. Davis, R. Westfall, A. Flint, M. Ikegami, Hongfang Wang, D. Grivet (2010)
Gene movement and genetic association with regional climate gradients in California valley oak (Quercus lobata Née) in the face of climate change
Molecular Ecology, 19
S. Morishita (1998)
On Classification and Regression
(2008)
Climate estimates and plantclimate relationships
P. Pearman, M. D’Amen, C. Graham, W. Thuiller, N. Zimmermann (2010)
Within‐taxon niche structure: niche conservatism, divergence and predicted effects of climate change
Ecography, 33
M. Lascoux, A. Palmé, R. Cheddadi, R. Latta (2004)
Impact of Ice Ages on the genetic structure of trees and shrubs.
Philosophical transactions of the Royal Society of London. Series B, Biological sciences, 359 1442
W.W. Oliver, R.A. Ryker, R.M. Burns, B.H. Honkala (1990)
Silvics of North America: 1. Conifers,
, vol 1
D. Rosauer, Renee Catullo, J. VanDerWal, A. Moussalli, C. Moritz (2015)
Lineage Range Estimation Method Reveals Fine-Scale Endemism Linked to Pleistocene Stability in Australian Rainforest Herpetofauna
PLoS ONE, 10
J. Norris, S. Jackson, J. Betancourt (2006)
Classification tree and minimum‐volume ellipsoid analyses of the distribution of ponderosa pine in the western USA
Journal of Biogeography, 33
S.E. Williams, L.P. Shoo, J.L. Isaac, A.A. Hoffmann, G. Langham (2008)
PLoS Biol.
, 6
O. Broennimann, M. Fitzpatrick, P. Pearman, B. Petitpierre, L. Pellissier, N. Yoccoz, W. Thuiller, M. Fortin, C. Randin, N. Zimmermann, C. Graham, A. Guisan (2012)
Measuring ecological niche overlap from occurrence and spatial environmental data
Global Ecology and Biogeography, 21
D.J. Shinneman, R.E. Means, K.M. Potter, V.D. Hipkins (2016)
PLoS One
, 11
F. Sorensen, N. Mandel, J. Aagaard (2001)
Role of selection versus historical isolation in racial differentiation of ponderosa pine in southern Oregon: an investigation of alternative hypotheses
Canadian Journal of Forest Research, 31
B. Oney, B. Reineking, G. O'neill, J. Kreyling (2013)
Intraspecific variation buffers projected climate change impacts on Pinus contorta
Ecology and Evolution, 3
Wells (1964)
Geographic variation in Ponderosa Pine.
R. Latta, Y. Linhart, David Fleck, M. Elliot (1998)
DIRECT AND INDIRECT ESTIMATES OF SEED VERSUS POLLEN MOVEMENT WITHIN A POPULATION OF PONDEROSA PINE
Evolution, 52
E. Rudloff, M. Lapp (1992)
Chemosystematic studies in the genus Pinus. VII. The leaf oil terpene composition of ponderosa pine, Pinus ponderosa
Botany, 70
Rehfeldt (1990)
Genetic differentiation among populations of Pinus ponderosa from the Upper Colorado
(1989)
Paternal inheritance of chloroplastDNA and maternal inheritance of mitochondrial DAN in loblolly pine
G. Rehfeldt (2018)
Adaptive Variation in Pinus Ponderosa From Intermountain Regions: II. Middle Columbia River System
R. Latta, J. Mitton (1999)
HISTORICAL SEPARATION AND PRESENT GENE FLOW THROUGH A ZONE OF SECONDARY CONTACT IN PONDEROSA PINE
Evolution, 53
Keller (2004)
Stress and respiration traits differ among four geographically distinct Pinus ponderosa seed sources. Thermochim.
(2016)
ecospat: spatial ecology miscellaneous methods
D.B. Neale, R.R. Sederoff (1989)
Theor. Appl. Genet.
, 77
F. Valladares, S. Matesanz, F. Guilhaumon, M. Araújo, L. Balaguer, M. Benito-Garzón, W. Cornwell, E. Gianoli, M. Kleunen, D. Naya, A. Nicotra, Hendrik Poorter, M. Zavala (2014)
The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change.
Ecology letters, 17 11
Oliver (1990)
Ponderosa pine.
M. Davis (1989)
Lags in vegetation response to greenhouse warming
Climatic Change, 15
R. Dumroese, Mary Williams, John Stanturf, J. Clair (2015)
Considerations for restoring temperate forests of tomorrow: forest restoration, assisted migration, and bioengineering
New Forests, 46
Andy Liaw, M. Wiener (2007)
Classification and Regression by randomForest

Publisher: Oxford University Press
Copyright: Copyright © 2022 Society of Systematic Biologists
ISSN: 1063-5157
eISSN: 1076-836X
DOI: 10.1093/sysbio/syy017
Publisher site: See Article on Publisher Site

Abstract

Abstract Unique responses to climate change can occur across intraspecific levels, resulting in individualistic adaptation or movement patterns among populations within a given species. Thus, the need to model potential responses among genetically distinct populations within a species is increasingly recognized. However, predictive models of future distributions are regularly fit at the species level, often because intraspecific variation is unknown or is identified only within limited sample locations. In this study, we considered the role of intraspecific variation to shape the geographic distribution of ponderosa pine (Pinus ponderosa), an ecologically and economically important tree species in North America. Morphological and genetic variation across the distribution of ponderosa pine suggest the need to model intraspecific populations: the two varieties (var. ponderosa and var. scopulorum) and several haplotype groups within each variety have been shown to occupy unique climatic niches, suggesting populations have distinct evolutionary lineages adapted to different environmental conditions. We utilized a recently available, geographically widespread dataset of intraspecific variation (haplotypes) for ponderosa pine and a recently devised lineage distance modeling approach to derive additional, likely intraspecific occurrence locations. We confirmed the relative uniqueness of each haplotype-climate relationship using a niche-overlap analysis, and developed ecological niche models (ENMs) to project the distribution for two varieties and eight haplotypes under future climate forecasts. Future projections of haplotype niche distributions generally revealed greater potential range loss than predicted for the varieties. This difference may reflect intraspecific responses of distinct evolutionary lineages. However, directional trends are generally consistent across intraspecific levels, and include a loss of distributional area and an upward shift in elevation. Our results demonstrate the utility in modeling intraspecific response to changing climate and they inform management and conservation strategies, by identifying haplotypes and geographic areas that may be most at risk, or most secure, under projected climate change. Ecological niche modeling, Pinus ponderosa, intraspecific variation, haplotype, range shift, conservation, climate-plant relationships A species’ response to environmental change is not always uniform across its distribution as genotypic and phenotypic differences within a species can interact with local environmental conditions and result in local adaptation (Valladares et al. 2014). Even though most biogeographical modeling studies treat a species as a single unit, its response to climate change will typically vary across space, either through plastic gene expression, adaptation, extirpation, or migration (Sork et al. 2010; Oney et al. 2013; Valladares et al. 2014). Therefore, population and phylogeographic structure should not be ignored when modeling species distributions, inferring niche structure, or predicting how climate will affect species distributions (Pfenninger et al. 2007; Pearman et al. 2010). Not including phylogeographic structure may lead to underrepresentation of certain lineages in species models (Pearman et al. 2010), especially when they are clearly defined and occupy distinct niche spaces (Hällfors et al. 2016). Distinct lineages defined by genetic variation are the units of evolutionary change and thus crucial for providing insight into the ability of a species to adapt to environmental change. Predicting intraspecific geographic response to environmental factors is becoming increasingly feasible as research reveals genetic variation across the distribution of individual species (e.g., Pearman et al. 2010; Espíndola et al. 2012; D’Amen et al. 2013; Yannic et al. 2014; Hällfors et al. 2016; Prasad and Potter 2017). Modeling the ecological niches of populations or lineages may be especially important for species with spatial variation in genotypic or phenotypic traits that suggest they are in the process of differentiating into two or more species (Valladares et al. 2014). Ponderosa pine (Pinus ponderosa) is a widespread, well-studied, and economically important species in Western North America, and its substantial morphological variation (Wells 1964) suggests it may be in early stages of differentiation into multiple species (Wang 1977, Jaramillo-Correa et al. 2009). Two varieties within the species are recognized—the Pacific variety, P. ponderosa var. ponderosa Laws., has open, plume-like foliage, and a low proportion of two-needle fascicles as opposed to the Rocky Mountain variety, P. ponderosa var. scopulorum Engelm., which has compact and brush-like foliage and moderate to high proportions of two-needle fascicles. The separation of the two varieties is also established by leaf oil terpene composition (Rudloff and Lapp 1992), mitochondrial DNA (Johansen and Latta 2003; Potter et al. 2013), highly polymorphic nuclear microsatellite and isozyme markers (Potter et al. 2015), and allozymes (Niebling and Conkle 1990). The two varieties have similar temperature limitations but occur within different precipitation regimes, with var. ponderosa locations generally corresponding to higher winter moisture and var. scopulorum to higher summer moisture availability (Norris et al. 2006; Shinneman et al. 2016). It is thought that the varieties diverged at least before the last glacial maximum, and more likely around 250,000 years ago (Lascoux et al. 2004), with Pleistocene glacial and interglacial cycles leading to further genetic diversification through isolation in localized refugia, followed by subsequent migration, hybridization, and introgression (Shinneman et al. 2016). Within each variety, several geographically structured haplotypes have also been identified using mitochondrial DNA (Johansen and Latta 2003; Potter et al. 2013) and plastome sequences (Wofford et al. 2014), and genetic clusters have been identified from nuclear DNA microsatellite markers (Potter et al. 2015). Geographic structure identified from mtDNA suggests each haplotype resulted from unique patterns of migration, isolation, and local adaptation (Potter et al. 2013). Mitochondrial DNA is maternally inherited and dispersed via seeds with limited dispersal compared with pollen-dispersed chloroplast DNA, and geographic structure in genetic differentiation is retained longer (Neale and Sederoff 1989; Petit et al. 1993). Thus, the haplotypes may represent evolutionarily distinct units capable of responding differently to climate change as a result of unique adaptations to refugial climate conditions and long-term isolation during glacial periods (Potter et al. 2013, Shinneman et al. 2016). Given unique intraspecific lineages and the potential for distinctive climate niches, the need to model distributions of ponderosa pine at the subspecific level has been recognized (Norris et al. 2006; Rehfeldt et al. 2014, Shinneman et al. 2016). Shinneman et al. (2016) modeled the contemporary and historical climate niche distributions of the two varieties, as well as 10 geographically structured mtDNA haplotypes identified by Potter et al. (2013) to further infer the evolutionary and phylogeographic history of these lineages. Recent modeling research has also suggested that ponderosa pine may be threatened by changing climate, but that the response may differ for the two varieties (Rehfeldt et al. 2014). However, the potential for future climate change to shape the distribution of ponderosa pine populations, as defined by haplotypic diversity, has not been explored, and has only been studied for a handful of other species relative to intraspecific variation (Pearman et al. 2010; Bálint et al. 2011; Bentio Garzón et al. 2011; Oney et al. 2013; Prasad and Potter 2017). Given past climate influences on ponderosa pine evolutionary trajectories and genetic lineages, we sought to better understand potential future distributions of the species by specifically considering responses to expected climate change of the two well-known and broadly distributed varieties and the haplotypes associated with them. Specifically, we were interested in determining how haplotype-level ecological niche models (ENMs) might differ from variety-level ENMs in terms of predicting recent and future geographic ranges in relation to climate, and exploring the advantages and challenges of estimating ENMs below the traditional level of species. To accomplish this objective, we employed a three-tiered approach. First, to address potentially inadequate geographic coverage of samples across populations, we utilized a lineage distance modeling approach (Rosauer et al. 2015) to derive likely haplotype occurrence locations (presences) by leveraging distribution information for the two varieties and known haplotype locations. Second, both actual and derived presence-absence datasets were employed in an ensemble of ENMs to project recent (1961–1990) and future (2060, 2090) distributions for ponderosa pine across three intraspecific levels: 1) for each individual haplotype; 2) for all haplotypes within a variety combined; and 3) for each corresponding variety. Third, to quantify the uniqueness of each estimated haplotype niche space relative to all others, we applied a post hoc, niche-overlap analysis. The combination of these methodological approaches allowed us to provide relatively robust projections of future distributions of ponderosa pine across intraspecific levels, and permitted a more nuanced interpretation of potentially unique population-level responses to climate change that may better inform conservation efforts for the species. Methods Haplotype Occurrence Data We examined eight of the ten haplotypes identified in Potter et al. (2013). We did not examine haplotypes 9 and 10, each limited to a small area of California, because of restricted and localized populations less amenable to niche modeling. Haplotypes were identified from the mtDNA nad1 second intron minisatellite region of foliage samples collected from 3113 trees representing 104 populations of ponderosa pine across its range (Potter et al. 2013). At least 20 individual trees were sampled in each population (except two), with most populations encompassing at least 30 sampled trees. Sampled trees were at least 100 m apart to encompass the entire range of the population’s genetic composition. Haplotypes 1, 5, and 8 are within the distribution of var. ponderosa and group together in phylograms of genetic distance (Potter et al. 2013). Haplotypes 2, 3, 4, 6, and 7 are within the distribution of var. scopulorum and are phylogenetically distinct from the haplotypes in var. ponderosa (Potter et al. 2013). Duplicate occurrences of the same haplotype that fell within the same geographic grid cell (800 m |$\times$| 800 m) were removed from the dataset. Climate Data Climate variables derived using thin-plate splines to interpolate 1961–1990 monthly temperature and precipitation data at 800 m resolution were obtained from the Research on Forest Climate Change website (charcoal.cnre.vt.edu/climate; Rehfeldt 2006; Crookston and Rehfeldt 2008), as well as an ensemble of 17 global climate models (Supplementary Appendix S1 available on Dryad at http://dx.doi.org/10.5061/dryad.jj264) for the RCP60 (medium-high greenhouse gas concentration pathway) scenario for 2060 and 2090. These climate layers were restricted to the western United States. Additional climate variables were derived from these climate variables (Table 1). We removed one of any two variables that was highly correlated (Pearson correlation |$>$| 0.7), keeping the more biologically relevant variable of a correlated pair (Rehfeldt et al. 2014; Shinneman et al. 2016). This resulted in seven variables used to fit the final models (except for random forest models of the varieties—see below): growing season precipitation, summer-spring precipitation balance, degree days below 0|$^{\circ}$|C, mean temperature of the warmest month, summer-winter temperature differential, annual moisture index, and precipitation ratio. Variables were log or square root transformed if needed to normalize data. Table 1. Climate variables considered for niche modeling of ponderosa pine varieties and haplotypes Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| MAP |$=$| mean annual precipitation; GSP |$=$| growing season precipitation; MTCM |$=$| mean temperature of the coldest month; MTWM |$=$| mean temperature of the warmest month. |$^{\mathrm{a}}$|All variables were used in an initial random forest model for variety distributions, and seven variables with Pearson correlations |$<$|0.7 were used in final niche models. |$^{\mathrm{b}}$|Derived variables are indicated. Open in new tab Table 1. Climate variables considered for niche modeling of ponderosa pine varieties and haplotypes Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| Climate variables . MAP |$^{\mathrm{a}}$|GSP; Apr. through Sept. Summer precipitation; July|$+$|Aug. Summer precipitation balance; (July|$+$|Aug.|$+$|Sept.)/(Apr.|$+$|May|$+$|June) |$^{\mathrm{a}}$|Summer–spring precipitation balance; (July|$+$|Aug.)/(Apr.|$+$|May) Spring precipitation; Apr.|$+$|May Winter precipitation; Nov. through Feb. Julian date the sum of degree days |$>$| 5|$^{\circ}$|C reaches 100 |$^{\mathrm{a}}$|Degree-days |$< 0^{\circ}$|C (based on mean monthly temperature) Degree-days |$>$| 5|$^{\circ}$|C Julian date of the first freezing date in autumn Length of frost-free period (days) Julian date of the last freezing date of spring Degree-days |$>$| 5|$^{\circ}$|C accumulating within the frost-free period (GSDD5) Degree-days > 0|$^{\circ}$|C (based on mean min. monthly temperature) MTCM |$^{\mathrm{a}}$|MTWM |$^{\mathrm{a}}$|Summer-winter temperature differential; MTWM-MTCM (TDIFF)|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Annual moisture index; (⁠|$\surd $|GSDD5)/MAP|$^{\mathrm{b}}$| Summer moisture index; (⁠|$\surd $|GSDD5)/GSP|$^{\mathrm{b}}$| |$^{\mathrm{a}}$|Precipitation ratio; GSP/MAP|$^{\mathrm{b}}$| TDIFF |$\times$| log(GSP)|$^{\mathrm{b}}$| Summer heat:moisture index; MTWM/log(May|$+$|June precip.)|$^{\mathrm{b}}$| MAP |$=$| mean annual precipitation; GSP |$=$| growing season precipitation; MTCM |$=$| mean temperature of the coldest month; MTWM |$=$| mean temperature of the warmest month. |$^{\mathrm{a}}$|All variables were used in an initial random forest model for variety distributions, and seven variables with Pearson correlations |$<$|0.7 were used in final niche models. |$^{\mathrm{b}}$|Derived variables are indicated. Open in new tab Lineage Distance Modeling for Haplotypes Although predicting the niches of intraspecific populations is more feasible given the rise in research that uncovers genetic variation of species across their distributions, often there are not enough sampled occurrences of each population to define their climate niche space and fit robust ENMs, limiting the ability to hindcast or forecast their distributions (Rosauer et al. 2015). In addition, it may not be appropriate to fit an ENM directly to only known haplotype occurrences because these occurrences may be governed by other factors than environmental conditions, such as competition. In the case of ponderosa pine, there may be too few sampled occurrences among intraspecific populations to produce robust ENMs, and sampled occurrences alone most likely do not represent the entire climatic niche of each haplotype. Thus, to overcome low sample sizes of known haplotype occurrences, we employed a recently developed method, lineage distance modeling (Rosauer et al. 2015), that probabilistically estimates the likely distribution of intraspecific lineages, therefore, increasing the number of occurrence points to fit ENMs. This is accomplished by first fitting an initial ENM for a species, and then assigning pixels within the species distribution to haplotypes based on the distance to known haplotype occurrences, taking into account barriers to dispersal using a measured cost distance (Fig. 1). The cost-distance analysis is calculated as |$-$|log(habitat suitability) where habitat suitability is the probability of occurrence from the species ENM. The output is a probability map for each lineage, such that all estimated probabilities for a cell sum to the total probability of the species occupying the cell from the original ENM. For this analysis, ENMs were fit for the two varieties separately (as described below), rather than the entire species, because they have distinct climate niches (Norris et al. 2006). The resulting distribution maps for each haplotype are effectively defined by the niche of the variety and geographically constrained by theoretical dispersal limitations. Custom python code from Rosauer et al. (2015) was used in ArcMap 10.4.1 to develop the haplotype probability maps. Estimated probability maps for each haplotype were compared with known haplotype locations using the Brier Score (Brier 1950), a measure of the accuracy of probabilistic predictions to binary (presence/absence) observations calculated as the mean squared difference between the two. The lower the Brier Score, the better the predictions. Figure 1. Open in new tabDownload slide Workflow of lineage distance modeling and ENM shown for var. ponderosa and haplotypes 1, 5, and 8. First, an ENM is fit for the variety using RF (Step 1). Then, lineage distance modeling is used to produce distribution maps for each haplotype (Step 2). These distribution maps are sampled, and ENMs are fit for each haplotype using RF, MARS, and NPMR (Step 3). Lastly, ENMs are fit for the variety using RF, MARS, and NPMR (Step 4). An ensemble distribution is created for each haplotype and the variety by averaging the RF, MARS, and NPMR models. Figure 1. Open in new tabDownload slide Workflow of lineage distance modeling and ENM shown for var. ponderosa and haplotypes 1, 5, and 8. First, an ENM is fit for the variety using RF (Step 1). Then, lineage distance modeling is used to produce distribution maps for each haplotype (Step 2). These distribution maps are sampled, and ENMs are fit for each haplotype using RF, MARS, and NPMR (Step 3). Lastly, ENMs are fit for the variety using RF, MARS, and NPMR (Step 4). An ensemble distribution is created for each haplotype and the variety by averaging the RF, MARS, and NPMR models. Variety and Haplotype Niche Modeling To fit ENMs for each variety, 20,000 cells were subsampled from ponderosa pine distribution maps of the USDA Forest Service’s National Individual Tree Species Atlas (https://www.fs.fed.us/foresthealth/applied_sciences/mapping-reporting/remote-sensing/index.shtml) that were also occupied by vegetation types that contain ponderosa pine in the U.S. Geological Survey National Gap Analysis Program (https://gapanalysis.usgs.gov) map (Fig. 1). Because both of these mapped sources were derived from classifications of satellite imagery, we used agreement between the two datasets for presence locations to minimize errors of commission (similar to Shinneman et al. 2016). Each of the occupied cells was assigned to one of the two varieties based on its geographic position. In transition zones between the two varieties (southern California and central Montana), the closest known haplotype occurrence was used to assign cells to a variety. Lastly, known haplotype occurrences were added to their corresponding variety. In the case of multiple occurrences of a single haplotype at a given location, only one record per that haplotype per grid cell (800 m |$\times$| 800 m) was kept for analysis. An initial 20,000 absence points were randomly subsampled, and 500 additional absences were generated for each variety from the cells occupied by the other variety. To fit ENMs for each haplotype, original haplotype presences from Potter et al. (2013) were augmented with up to 1000 cells sampled from the distribution map derived from the lineage distance modeling method (Fig. 1). Sampling from the lineage distance map was weighted so that 75% of samples had probabilities between 0.9 and 1, and 25% were between 0.5 and 0.89. Although this weighting scheme is a conservative approach to estimating haplotype occurrences and minimizes potential overestimation in peripheral populations, the lineage distance modeling also generally added samples in peripheral locations relative to the actual haplotype locations (Supplementary Appendix S2 available on Dryad), potentially balancing out any bias against peripheral populations. Indeed, after testing two less conservative weighting approaches, we found no difference between the ability of each sampling scheme to predict known haplotype locations (i.e., Brier Scores were not significantly different). Moreover, we also explored the possibility that using raw probability of occurrence values from the lineage distance modeling (analogous to using abundance) would produce less-restrictive ENMs for haplotypes, but the projected distributions were similar (or even more restricted, see Supplementary Appendix S3 available on Dryad). We also found little difference in Brier Scores using raw probability of occurrence values compared with the threshold sampling procedure (Supplementary Appendix S3 available on Dryad). Absences for each haplotype consisted of original absences from Shinneman et al. (2016) and absences generated in the previous step for the variety to which it belongs. We used the absences of each variety rather than of the entire species because the delineation of haplotypes between the varieties is robust based on mtDNA and nuclear DNA analyses (Potter et al. 2013; Potter et al. 2015). For all haplotype and variety niche models, data were split into 10 random partitions of training (80%) and testing (20%) sets. To minimize uncertainty in model algorithm choice, three different ENM modeling approaches were used to fit each partition (i.e., resulting in 30 models per haplotype/variety): random forest (RF), multivariate adaptive regression spline (MARS), and non-parametric multiplicative regression (NPMR). RF and MARS models were fit and projected in R v 3.3.1 (R Core Team 2014) using the randomForest (Liaw and Wiener 2002) and earth (Milborrow 2016) packages, respectively. NPMR models were run in HyperNiche 2.28 (McCune and Mefford 2009). All 23-climate variables were input into the initial RF models built for the varieties that provided inputs for the lineage distance modeling to explore all climate variables and minimize under-prediction of variety distributions (RF randomly chooses predictor variables and eliminates those that are not improving model performance; therefore, removal of correlated variables a priori is not needed.) However, to facilitate compatibility of results among different modeling approaches, final variety and haplotype RF models were fit with the same seven uncorrelated climate variables used in the MARS and NPMR models (Table 1). Each individual haplotype and variety model was projected onto current and future climatic scenarios, and an ensemble model was created by first averaging the 10 model runs of ecological niche maps estimated from each modeling approach, and then averaging the three resulting maps per model approach for each climate scenario. We also combined the distributions of haplotypes within each variety (H1, H5, and H8 for var. ponderosa; H2, H3, H4, H6, and H7 for var. scopulorum) at each time period in order to observe how the distribution of the variety produced from ENM haplotypes compared with the distribution produced by an ENM fit with variety occurrences. These combined distribution maps represent the probability of at least one of the haplotypes occurring in a particular cell (i.e., by calculating |$1 - X$|⁠, where |$X$| is the probability that all haplotypes are absent). Model performance for each ENM (haplotypes and varieties) was tested using the area under the receiver operating curve (AUC), a metric used when comparing continuous outputs (ENM occurrence probabilities) to binary observations (presence/absence dataset sampled from the lineage distance method), ranging from 0 to 1 in which a value of 1 indicates perfect agreement and a value of 0.5 indicating model performance is no better than random. We also evaluated model performance by comparing the known haplotype locations to the predicted occurrence probabilities using the Brier Score. Niche-Overlap Analyses The amount of niche overlap and niche similarity between haplotypes was calculated following Broennimann et al. (2012), using the “PCA-env option”, in which a Kernel smoother is applied to densities of species occurrence in a gridded environmental space calibrated on the available environmental space. The D metric (Schoener 1970) measures the amount of niche overlap between two haplotypes in the gridded environmental space (D |$=$| 1: complete overlap; D |$=$| 0: no overlap): \[ D=1-0.5 \left(\sum_{ij} \left| {z1}_{ij}-{z2}_{ij} \right|\right) \] where |${z1}_{ij}$| and |${z2}_{ij}$| are the occupancies of entity 1 and 2, respectively, and |$i$| and |$j$| refer to the cell corresponding to the |$i$|th and |$j$|th bins of the environmental variables. The similarity test measures the amount of overlap of the niche of one haplotype with a random subset in the second haplotype’s range, and vice versa, for 100 iterations, and compares this distribution of overlap to the observed D metric. The similarity test can indicate if two niches are significantly similar, significantly dissimilar, or nonsignificant (indicating not enough statistical power to determine either/or). Niche space here was defined using the same presence points used in ENM modeling. Niche overlap and similarity analyses were performed using the ecospat package (Broennimann et al. 2016) in R. Results Lineage Distance Modeling for Haplotypes Distribution maps of each haplotype from the lineage distance modeling are generally in agreement with known locations of each haplotype (Brier Scores |$<$| 0.2; Table 2; Supplementary Appendix S2 available on Dryad). However, there were limitations to this approach for certain haplotypes. Specifically, because the distributions of H2, H4, and H7 are geographically restricted (Fig. 3 in Supplementary Appendix S2 available on Dryad), fewer than 1000 occurrence points from the lineage distance maps for these haplotypes were available for ENM construction (Table 2; H4: |$n =$| 102, H2: |$n =$| 296, and H7: |$n =$| 253). Furthermore, because the ENM for var. scopulorum predicted only low probabilities of occurrence (⁠|$<$| 0.50) for the small and isolated populations of H2 in southern California/Nevada and H7 in Nevada, in those specific areas only the original haplotype locations (from Potter et al. 2013) were used to construct ENMs. Figure 3. Open in new tabDownload slide Percent area change and mean elevation change from the recent period to 2060 and 2090 for the varieties (Var.), combined haplotypes (Haps.), and individual haplotypes. Change values are calculated among time periods using cells with probability of occurrence estimates |$>$| 0.50. Figure 3. Open in new tabDownload slide Percent area change and mean elevation change from the recent period to 2060 and 2090 for the varieties (Var.), combined haplotypes (Haps.), and individual haplotypes. Change values are calculated among time periods using cells with probability of occurrence estimates |$>$| 0.50. Table 2. Lineage distance modeling and ensemble ENM evaluation Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Notes: Brier scores represent the mean squared difference between the predicted probabilities of occurrence and the observed binary presence/absence of each haplotype for the lineage distance method and ensemble ENM. Mean AUC values across 10 model runs for all three models for each variety and haplotype indicate how well ENMs discriminates between presences and absences. Lastly, the number of presences for each variety and haplotype used to build the ENMs are listed, with the number of actual locations (from Potter et al. 2013) followed by estimated locations (using the linear distance method) in parentheses. Open in new tab Table 2. Lineage distance modeling and ensemble ENM evaluation Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Variety/ Haplotype . Lineage distance Brier score . ENM Brier score . ENM AUC . Number of presences . var. ponderosa — — 0.98 104(12,230) var. scopulorum — — 0.97 144(7770) H1 0.11 0.20 0.97 25(998) H2 0.05 0.03 0.99 18(296) H3 0.14 0.19 0.96 64(998) H4 0.03 0.02 0.98 14(102) H5 0.04 0.10 0.99 25(939) H6 0.11 0.17 0.95 52(620) H7 0.07 0.10 0.97 29(253) H8 0.06 0.14 0.97 26(998) Notes: Brier scores represent the mean squared difference between the predicted probabilities of occurrence and the observed binary presence/absence of each haplotype for the lineage distance method and ensemble ENM. Mean AUC values across 10 model runs for all three models for each variety and haplotype indicate how well ENMs discriminates between presences and absences. Lastly, the number of presences for each variety and haplotype used to build the ENMs are listed, with the number of actual locations (from Potter et al. 2013) followed by estimated locations (using the linear distance method) in parentheses. Open in new tab Variety and Haplotype Niche Models ENMs for varieties and haplotypes were able to discriminate well between occupied and unoccupied cells (AUC scores |$>$| 0.9; Table 2), and predicted the occurrence of known haplotype occurrences well (Brier Score |$< = 2.0$|⁠, Table 2). Variable importance differed for each variety and haplotype as well as among model type (Supplementary Appendix S4 available on Dryad). Var. ponderosa.— The ensemble ENM for the variety matches its known distribution, although it slightly over-predicts in portions of the Great Basin (in Nevada, Southern Idaho, and Eastern Oregon) and the Wasatch Range in Northern Utah (Fig. 2a). Projections of the ENM for the variety to 2060 and 2090 (Fig. 2b and c) suggest a range contraction will occur, with a projected 24% loss in total area within high occurrence probability areas (⁠|$>$| 0.5 occurrence probability) by 2090 (compared with the recent period) (Fig. 3), especially at lower elevations in central Oregon, northeastern Washington, and northeastern California. This trend is compensated slightly with potential range extension into higher elevations, especially in central Idaho and the southern Sierra Nevada of California. Indeed the mean elevation for cells with high probability (⁠|$>$| 0.5) of occurrence is projected to increase by 132 m on average by 2090 (compared with the recent period) (Fig. 3). Moreover, the mean probability of all occupied cells with |$>$| 0.5 probability in the recent period is generally projected to decrease by 0.02 and then by 0.53 in 2060 and 2090, respectively (Fig. 4). The combined distribution for the three haplotypes (H1, H5, and H8) within the variety is not as extensive compared with the variety model for the recent climate period (Fig. 2d), but similar distributional trends, such as loss in central Oregon, are projected through 2090 (Fig. 2e–f). However, for the combined haplotype distribution there is no loss in total area and only a 57 m increase in mean elevation for cells with high probability (⁠|$>$| 0.5) of occurrence (Fig. 3), although mean probability values for all occupied cells in the recent period with |$>$| 0.5 probability occurrence are projected to decline for both 2060 and 2090, by 0.17 and 0.53, respectively (Fig. 4). Individual haplotype models generally project similar trends to the variety (Supplementary Appendix S5 available on Dryad), as highlighted here for H1 (Fig. 2g–i), which is projected to have a decrease in area (41%), an increase in mean elevation (132 m), and a decrease in mean probability (by 0.63, of occupied cells in the recent period with |$>$| 0.5 probability) for 2090 (Figs. 3 and 4). Figure 2. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. ponderosa (a–c), combined haplotype (H1, H5, and H8) models (d–f), and a single haplotype (H1) model (g–i). Predicted distributions in all cases are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. ponderosa (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Figure 2. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. ponderosa (a–c), combined haplotype (H1, H5, and H8) models (d–f), and a single haplotype (H1) model (g–i). Predicted distributions in all cases are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. ponderosa (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Figure 4. Open in new tabDownload slide Mean occurrence probability of occupied cells (with |$>$| 0.5 occurrence probability) from ensemble model of recent distributions of ponderosa pine lineages, combined haplotypes, and individual haplotypes and the mean occurrence probability of those cells in 2060 and 2090. Figure 4. Open in new tabDownload slide Mean occurrence probability of occupied cells (with |$>$| 0.5 occurrence probability) from ensemble model of recent distributions of ponderosa pine lineages, combined haplotypes, and individual haplotypes and the mean occurrence probability of those cells in 2060 and 2090. Var. scopulorum.— The modeled niche distribution of var. scopulorum matches the known distribution of the variety, although there is some over-prediction of the distribution in Montana and under prediction in the Great Plains (Fig. 5a). In general, the ensemble ENM of the variety and the individual haplotypes (H3 and H6) could not accurately model the northern portion of the variety’s distribution. The easternmost populations of the variety along the Niobrara River in Nebraska are not predicted by the ensemble ENMs most likely because these populations exist in microclimates along the river that are not available in surrounding areas, making it difficult for coarse scale climate data (800 m x 800 m resolution) to capture their presence. Figure 5. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. scopulorum (a–c), combined haplotype (H2, H3, H4, H6, and H7) models (d–f), and a single haplotype (H3) model (g–i). Predicted niche distributions are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. scopulorum (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Figure 5. Open in new tabDownload slide Predicted niche distributions of lineages of ponderosa pine under recent (1961–1990) and projected future (2060, 2090) climates for var. scopulorum (a–c), combined haplotype (H2, H3, H4, H6, and H7) models (d–f), and a single haplotype (H3) model (g–i). Predicted niche distributions are from an ensemble of three ENMs, and future climates are ensemble projections (see Methods section). Black outline indicates approximate recent distribution of var. scopulorum (Little 1971), and overlain known haplotype occurrence locations are from Potter et al. (2013). Similar to var. ponderosa, ENM projections for var. scopulorum in the future suggest a range contraction, with a 59% area loss and a 670 m increase in mean elevation by 2090 compared with the recent time period for areas with |$>$| 0.5 occurrence probability (Fig. 3). There is some predicted expansion northward into central Utah and western Colorado, and some expansion into southwest Montana even though there is predicted loss for most Montana populations. Occupied cells with probabilities |$>$| 0.5 in the recent period decline by 0.17 in 2060 and by 0.57 in 2090 (Fig. 4). The combined predicted niche distributions of the individual haplotypes within the variety (H2, H3, H4, H6, and H7) are more delimited than the variety (similar to the pattern observed with var. ponderosa haplotypes), but they have a similar pattern of loss and upward range shift as the variety, with a 23% loss in area and a mean elevational increase of 363 m for high probability (⁠|$>$| 0.5) occurrences by 2090 (Fig. 3). The probability of the combined haplotypes occurring in predicted occupied sites in the recent time period with |$>$| 0.5 probability in the recent time period decreases by 0.61 in 2090 (Fig. 4). Individual haplotypes within var. scopulorum follow a similar trend (Supplementary Appendix S5 available on Dryad), as highlighted by H3 (Fig. 5g–i), which is projected to have a decrease in area (40%), an increase in mean elevation (422 m), and a decrease in mean probability (by 0.63, of occupied cells in the recent period with |$>$| 0.5 probability) for 2090 (Figs. 3 and 4). Niche Overlap Haplotypes of ponderosa pine occupy distinct niche spaces in western North America, as no two haplotypes had statistically similar niche spaces (Similarity Test |$P \,{<}\,0.01$|⁠; Supplementary Appendix S6 available on Dryad). Thus, even though the niche spaces of some haplotypes overlap (e.g., H3 and H7), they are not significantly similar (Fig. 6 and Supplementary Appendix S6 available on Dryad). Haplotypes with larger distributional ranges tend to occupy more niche space (e.g., H1 vs. H4; Fig. 6; Supplementary Appendix S6 available on Dryad). Those haplotypes within var. ponderosa occupy environmental space that contains higher winter precipitation than haplotypes within var. scopulorum, in agreement with previous analyses (Norris et al. 2006; Shinneman et al. 2016; Fig. 6). Within var. ponderosa, H1 and H8 have the most overlap (D |$=$| 0.497), whereas H8 and H5 have very little overlap (D |$=$| 0.03). Within var. scopulorum, the most overlap is between H3 and H7 (D |$=$| 0.324) and H3 and H2 (D |$=$| 0.244) whereas H4 is the most isolated (D |$<$| 0.05 with all haplotypes within the variety). The niche space of H4, however, overlaps with the niche space of H1 (D |$=$| 0.126) and H8 (D |$=$| 0.260) of var. ponderosa (Fig. 6). Otherwise, haplotypes from var. ponderosa have very little niche overlap with haplotypes within var. scopulorum. H2, which has populations in New Mexico and southern California, has been placed within var. scopulorum (Potter et al. 2013) and the niche-overlap results support this (indeed, H2 overlaps with H3 [D |$=$| 0.244] more than other haplotypes). However, the low niche overlap between H2 and geographically proximate haplotypes in California associated with var. ponderosa, H1 (D |$=$| 0.023) and H5 (D |$=$| 0.014), could also be influenced by the lack of H2 presences predicted for southern California by the lineage distance modeling. Figure 6. Open in new tabDownload slide a) Niche of each variety and haplotype in climatic space defined by the two first axes of the PCA of the seven climate variables used in ENMs for western United States. The amount of variation explained by PC Axis 1 and 2 is 38.54% and 28.84%, respectively. Gray shading represents the density of occurrences of each variety and the lines indicate the outline of the niche space for each variety and haplotype. Inset is the contribution of the seven climate variables on the two axes of the PCA. b) The amount of niche overlap (D statistic) and |$P$|-values for niche similarity between each set of haplotypes (Supplementary Appendix S6 available on Dryad). Figure 6. Open in new tabDownload slide a) Niche of each variety and haplotype in climatic space defined by the two first axes of the PCA of the seven climate variables used in ENMs for western United States. The amount of variation explained by PC Axis 1 and 2 is 38.54% and 28.84%, respectively. Gray shading represents the density of occurrences of each variety and the lines indicate the outline of the niche space for each variety and haplotype. Inset is the contribution of the seven climate variables on the two axes of the PCA. b) The amount of niche overlap (D statistic) and |$P$|-values for niche similarity between each set of haplotypes (Supplementary Appendix S6 available on Dryad). Discussion Comparisons Among Projected ENM Distributions for Varieties Versus Haplotypes Phylogeographic structures defined by genetic variation may represent clearly defined and distinct niche spaces among lineages within a species that reflect the evolutionary capacity of a species to adapt to environmental change (Prasad and Potter 2017). Thus, modeling the ecological niches of intraspecific lineages may project distributions of each lineage that differ in important ways from projections for the species (Hällfors et al. 2016). For instance, using the fundamental niche of simulated conceptual species and subpopulations, Valladares et al. (2014) found more restricted distributions when modeling distinct subpopulations compared with modeling the entire species. In contrast, Pearman et al. (2010), Oney et al. (2013), and Valladares et al. (2014) found the opposite when modeling the realized niches of real species and populations within them. Valladares et al. (2014) attributed such differences to modeling the fundamental versus realized niches. However, our work suggests an alternative explanation is in order, as our realized niches of the haplotypes were generally more limiting than those of the varieties. Estimates derived from the collective representation of the haplotype probabilities for the recent period are also generally more geographically restricted and have higher average occurrence probabilities than those projected for the varieties (Figs. 2, 4, and 5), and these more restricted ranges are further amplified for most individual haplotypes under future climate change compared with the varieties (Supplementary Appendix S5 available on Dryad). Although this could be partly caused by the influence of haplotype presence locations that restrict estimates of the peripheral climate niche, we suggest that another plausible explanation is that the haplotype populations correspond to generally well defined and often minimally overlapping climate niches (Fig. 6) that are uniquely projected to decrease in the future. In contrast, models for the varieties capture a broader niche representing more potential combinations of climate conditions, a larger portion of which may persist into the future. It has been suggested that by not including intraspecific variation in ENMs, fitness-environment curves may be smoothed due to averaging across the broader niche, leading to more pessimistic forecasts of future distributions (Valladares et al. 2014). However, our variety models did not result in more unfavorable distributions (i.e., better delineated with higher probabilities) relative to the haplotype models. Regardless of intraspecific emphasis, by examining relative changes across models (based on probability estimates |$>$| 0.5), we were able to reveal potential trends that generally hold for both the variety and most corresponding haplotype models under climate change. For the Pacific variety (var. ponderosa), there is a general loss of area projected under climate change (by ~ 25%) and a modest increase in elevation (⁠|$<$| 200 m). Although individual haplotype responses were more variable, they generally supported trends of decreasing range size and increasing elevation projected by the variety model under future climates; however, the combined haplotype distributions did not (Figs. 2 and 3). Projected trends for the Rocky Mountain variety (var. scopulorum) were similar but more extreme (⁠|$>$| 50% decrease in area, and |$>$| 600 m elevation gain) than the Pacific variety, and were also generally supported by the more highly variable individual haplotype models, and to a lesser degree by the combined haplotype models (Figs. 3 and 5). These trends suggest that the individual haplotype models effectively captured inherent sub-specific variability in response, but that the combined haplotypes were less likely to decline overall because their ranges were already more restricted for the recent period (relative to ranges for the varieties), or because they were buoyed by a few favorable individual haplotype responses to climate change (e.g., projected distributions for H5, see Discussion below). We expected some predicted distribution loss in the future, especially given that extrapolation of ENMs fit with present-day climates and projected to novel climates tend to under-predict distributions (Maguire et al. 2016). This factor could be further enhanced when ENMs are fit with a narrower subset of occurrences, restricting potential future distributions. However, climate novelty is likely minimal in our models (⁠|$<$| 0.0003 following methods in Maguire et al. 2016), and the varying degrees of predicted area loss for haplotypes and varieties suggest potentially true signals in their projected responses to climate change. Indeed, the potentially unique responses of individual haplotypes to future climate within each ponderosa pine variety can be instructive, especially in terms of patterns of area loss and elevational shifts. Given likely long-term isolation and local evolution of each of the haplotypes within glacial periods, each may be expected to exhibit different responses to climate change (Potter et al. 2013; Shinneman et al. 2016). For example, within var. ponderosa, H5 does not lose a significant amount of area in the future under climate change (though probability values in estimated recent locations are projected to decline; Fig. 4), and it is the only haplotype projected to maintain a relatively stable mean elevation over time (Fig. 3). This may be due to ensemble climate projections that suggest relatively moderate warming and wet growing season conditions in the predominantly coastal and eastern Sierra Nevada range of this haplotype. Rehfeldt et al. (2014) also predicted relative geographic stability for var. ponderosa under future climate scenarios in this portion of its range. In contrast, within var. scopulorum, H6 is projected to lose most of its current range and shift significantly upwards (⁠|$>$| 1200 m) in mean elevation, most likely due to projections of hotter and drier growing seasons in the central and northeastern Rocky Mountains. The dramatic loss of H6 and var. scopulorum in the north (the lower elevations of the Black Hills and the hills and tablelands of south-eastern and central Montana) is in agreement with previous studies that modeled the niche distribution of the variety (Rehfeldt et al. 2014). The variety already occupies higher elevational areas in much of this part of its range, and there are large distances to the nearest analogous climate in the future, effectively restricting both elevational and latitudinal range shifts (Rehfeldt et al. 2014, Roberts and Hamann 2016). Haplotype 6 may also be vulnerable because it occupies a unique climate niche space (e.g., occupies higher values of growing season precipitation, precipitation ratio, and growing degree days below 0|$^{\circ}$|C than other haplotypes) that is on the edge of the total niche space for var. scopulorum (Fig. 6 and Supplementary Appendix S6 available on Dryad). Relevance to Adaptation and Phenotypic Plasticity To successfully adapt to climate change, plant species must either migrate to track their ecological niches or increase their tolerance to new climate conditions in situ. This latter strategy requires either expression of phenotypic plasticity, or selection on existing phenotypes defined by existing genetic material or on novel phenotypes resulting from favorable mutations (Williams et al. 2008; Chevin et al. 2010). Our haplotype ENMs operated under an assumption of adaptation to specific, realized climate niches; they did not directly incorporate plasticity or genetic adaptation under changing climate, and thus did not capture the ability of each haplotype to expand into or occupy other parts of the species’ fundamental niche. Moreover, we caution that the haplotypes, which are classified based on identification of neutral genetic markers, are not necessarily effective surrogates for uniquely adapted populations (e.g., Holderegger et al. 2006). However, although these mtDNA defined haplotypes do not necessarily have adaptive significance, they correspond to populations that have been influenced by long-term biogeographical processes. Previous studies of ponderosa pine have shown the existence of adaptive phenotypic variation within populations (e.g. Conkle 1973; Rehfeldt 1986a, 1986b, 1990; Namkoong and Conkle 1976; Sorensen and Weber 1994; Sorensen et al. 2001; Keller et al. 2004), and populations with greater phenotypic plasticity may be more tolerant of changing climates (Chevin et al. 2010). There are areas within the distribution of ponderosa pine with high genetic variation (based on nuclear DNA, see Potter et al. 2015), and greater genetic variation can allow populations to adapt to climate change via natural selection (Sork et al. 2010). In addition, high gene flow has been detected within each ponderosa pine variety, occurring across complex structures of genetic variation (Potter et al. 2015). Gene flow between even the most genetically structured populations can reduce population differentiation and relax migrational constraints (Potter et al. 2015). Given the potential for plasticity or genetically adaptive responses of ponderosa pine to changing climate, as well as shifts in biotic interactions, population-level responses may be best captured by developing ENMs for the varieties. Indeed, when reduced to two dimensions via the niche-overlap analysis, the collective ponderosa pine haplotype climate niche space does not cover the entire niche space of each variety (Fig. 6), potentially suggesting that the modeled niches of the varieties could indirectly represent more plasticity or adaptive potential for the species than the individual haplotypes combined. The ENMs of varieties may also be indirectly capturing intraspecific variation, and therefore, the plasticity of the variety, reflected in the projected broader geographic distributions and higher probabilities under climate change compared with most haplotypes. This highlights the complementary value of modeling at both the variety and haplotype levels, in relation to phenotypic plasticity and adaptation respectively, and providing a better understanding of how variability in intraspecific response might affect the future niche distribution of the species as a whole (Pauls et al. 2013; Marcer et al. 2016). Dispersal, Migration, and Competition Even though optimal climate for ponderosa pine is generally predicted to shift to higher elevations in the future, existing populations may not fill these areas because of dispersal limitations and migration lags (Davis 1989). Thus, rates of migration on leading edges may be slower than rates of habitat loss on trailing edges of existing populations (Campbell and Shinneman 2017). Yet, although ponderosa pine has limited optimal seed dispersal distance for regeneration (15–30 m; Oliver and Ryker 1990; Boyden et al. 2005; Latta et al. 1998), it is capable of low frequency, long-distance dispersal of 75–200 km (Lesser and Jackson 2012, 2013). Indeed, this likely explains its rapid expansion and northward migration from refugial areas during the Last Glacial Maximum, as supported by fossil evidence (Norris et al. 2016) and by a pattern of decreasing unique and rare alleles with increasing latitude found in contemporary populations (Potter et al. 2015). These paleogeographic range shifts may suggest potential for a similar response to future climate change. However, requirements of future migrations may be different from the past due to novel climate conditions, different refugia locations, and new migrational constraints imposed by human activity (Norris et al. 2016; Roberts and Hamann 2016). In addition, because realized climate niches (from which the ENMs were built) may be truncated due to biotic interactions; shifting interspecific competition under climate change could also affect future distributions (Wisz et al. 2013). Potential Genetic Variation and Diversity Under Future Climate Genetic variation and diversity are likely to be affected differently for ponderosa pine populations under future and potentially novel climate conditions compared with paleoclimates. For instance, based on fossil and genetic evidence (Potter et al. 2013), the contact zone between the two varieties in central Montana is relatively recent, established within the previous 1500 years, and likely controlled by a geographic shift from summer- to winter-dominant precipitation (Norris et al. 2006). Chloroplast and nuclear DNA patterns suggest introgression from west to east in this region (Latta and Mitton 1999; Potter et al. 2015), but there was no exchange of mtDNA (Latta and Mitton 1999; Johansen and Latta 2003; Potter et al. 2013). Potential genetic exchange between haplotypes and varieties in this location could be further restricted under future climate, based on the lack of overlap projected in our models for either the varieties or haplotypes. In addition, some haplotypes could become increasingly isolated under warmer and drier climates (e.g., H4 and H7). Thus, some genetic diversity may be lost, especially if isolated lineages have less suitable habitat and experience even greater geographic isolation under future climates. This underscores the need to identify potential future refugia across intraspecific levels that may offer the best chance for maintenance of genetic diversity and survival of the species under climate change (e.g., Keppel et al. 2012). Conclusions Our results highlight the potential value of modeling intraspecific response to climate change and, more specifically, are informative for future management and conservation of ponderosa pine, an economically, and ecologically important tree species. Ponderosa pine haplotype groups may represent evolutionary lineages adapted to different environmental conditions. Thus, by identifying lineages most vulnerable to climate change and detecting potential climatically suitable areas in the future, our results suggest at least three important considerations for successful conservation. First, it is significant that each haplotype occupies a statistically different environmental space within the range of ponderosa pine, as this could be considered evidence for differential adaptation to diverse environmental conditions and, thus, that the haplotypes may be evolutionary units worthy of distinct conservation attention. However, there is need for further study on adaptive genetic variation in ponderosa pine that could help to identify evolutionarily significant units for conservation, including defining genetic groups differently (i.e., based on nuclear DNA variation; Potter et al. 2015). Second, the results indicate that some haplotype groups, specifically H8 in var. ponderosa and H6 in var. scopulorum, could lose well more than half of their predicted distribution before 2090. These haplotypes may require special efforts to conserve their existing genetic variation, using both ex situ and in situ strategies (Potter et al. 2017). Peripheral populations of isolated lineages, such as H4 and H7, may also be at high risk under climate change, and due to their lower genetic diversity but greater genetic differentiation compared with populations in the core range of each ponderosa pine variety (Potter et al. 2015), may be ideal candidates for conservation action. Finally, the results of this study could help guide assisted species and population migration efforts (Dumroese et al. 2015; Aitken and Bemmels 2016), by identifying locations in which vulnerable ponderosa pine lineages may be best adapted in the future. Supplementary Material Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.jj264. Funding This work was supported in part by Cost Shares Agreement [14-CS-11330110-042 and 15-CS-11330110-067 to K.M.P.] between the Southern Research Station of the USDA Forest Service and North Carolina State University. Funding for developing the climate niche models and niche-overlap analysis presented in this article was provided by an Interagency Agreement [L12PG00219 to D.J.S.] between the U.S. Bureau of Land Management and the U.S. Geological Survey. Acknowledgments Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government. We would like to thank Dan Rosauer, Nicholas Matzke, and Brittany Barker for their very helpful reviews of this article. We would also like to acknowledge Robert Means, whose efforts helped to initiate this research, and who passed away before the work was completed. References Aitken S.N, Bemmels J.B. 2016 . Time to get moving: assisted gene flow of forest trees. Evol. Appl. 9 : 271 – 290 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Bálint M., Domisch S., Engelhardt C.H.M., Haase P., Lehrian S., Sauer J., Theissinger K., Pauls S.U., Nowak C. 2011 . Cryptic biodiversity loss linked to global climate change. Nat. Clim. Change 1 : 313 – 318 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Benito Garzón M., Alía R., Robson T.M., Zavala M.A. 2011 . Intra-specific variability and plasticity influence potential tree species distributions under climate change: intra-specific variability and plasticity. Global Ecol. Biogeogr. 20 : 766 – 778 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Boyden S., Binkley D., Shepperd W. 2005 . Spatial and temporal patterns in structure, regeneration, and mortality of an old-growth ponderosa pine forest in the Colorado Front Range. Forest Ecol. Manag. 219 : 43 – 55 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Brier G.W. 1950 . Verification of forecasts expressed in terms of probability. Mon. Weather Rev. 78 : 1 – 3 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Broennimann O., Fitzpatrick M.C., Pearman P.B., Petitpierre B., Pellissier L., Yoccoz N.G., Thuiller W., Fortin, M.-J., Randin C., Zimmermann N.E., Graham C.H., Guisan A. 2012 . Measuring ecological niche overlap from occurrence and spatial environmental data: Measuring niche overlap. Global Ecol. Biogeogr. 21 : 481 – 497 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Broennimann O., Di Cola V., Guisan A. 2016 . ecospat: spatial ecology miscellaneous methods. R package version 2.1.1. Available from: URL https://CRAN.R-project.org/package=ecospat. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Campbell J.L., Shinneman D.J. 2017 . Potential influence of wildfire in modulating climate-induced forest redistribution in a central Rocky Mountain landscape. Ecol. Process. 6 : 7 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Chevin L.-M., Lande R., Mace G.M. 2010 . Adaptation, plasticity, and extinction in a changing environment: towards a predictive theory. PLoS Biol. 8 : e1000357 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Conkle M.T. 1973 . Growth data for 29 years from the California elevational transect study of ponderosa pine. Forest Sci. 19 : 31 – 39 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crookston N.L., Rehfeldt G.E. 2008 . Climate estimates and plant-climate relationships. Moscow (Idaho) : USDA Forest Service, Rocky Mountain Research Station. Available from: URL http://forest.moscowfsl.wsu.edu/climate/. Google Scholar D’Amen M., Zimmermann N.E., Pearman P.B. 2013 . Conservation of phylogeographic lineages under climate change: African mammals and global warming. Global Ecol. Biogeogr. 22 : 93 – 104 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Davis M.B. 1989 . Lags in vegetation response to greenhouse warming. Climatic Change 15 : 75 – 82 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Dumroese R.K., Williams M.I., Stanturf J.A., St. Clair J.B. 2015 . Considerations for restoring temperate forests of tomorrow: forest restoration, assisted migration, and bioengineering. New Forests 46 : 947 – 964 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Espíndola A., Pellissier L., Maiorano L., Hordijk W., Guisan A., Alvarez N. 2012 . Predicting present and future intra-specific genetic structure through niche hindcasting across 24 millennia: hindcasting-based phylogeography. Ecol. Lett. 15 : 649 – 657 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hällfors M.H., Liao J., Dzurisin J., Grundel R., Hyvärinen M., Towle K., Wu G.C., Hellmann J.J. 2016 . Addressing potential local adaptation in species distribution models: implications for conservation under climate change. Ecol. Appl. 26 : 1154 – 1169 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Holderegger R., Kamm U., Gugerli F. 2006 . Adaptive vs. neutral genetic diversity: implications for landscape genetics. Landsc. Ecol. 21 : 797 – 807 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Jaramillo-Correa J.P., Beaulieu J., Khasa D.P., Bousquet J. 2009 . Inferring the past from the present phylogeographic structure of North American forest trees: seeing the forest for the genes. Can. J. Forest Res. 39 : 286 – 307 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Johansen A.D., Latta R.G. 2003 . Mitochondrial haplotype distribution, seed dispersal and patterns of postglacial expansion of ponderosa pine. Mol. Ecol. 12 : 293 – 298 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Keller E.A., Anekonda T.S., Smith B.N., Hansen L.D., Clair J.B.S., Criddle R.S. 2004 . Stress and respiration traits differ among four geographically distinct Pinus ponderosa seed sources. Thermochim. Acta 422 : 69 – 74 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Keppel G., Van Niel K.P., Wardell-Johnson G.W., Yates C.J., Byrne M. Mucina L., Schut A.G.T., Hopper S.D., Franklin S.E. 2012 . Refugia: identifying and understanding safe havens for biodiversity under climate change. Global Ecol. Biogeogr. 21 : 393 – 404 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lascoux M., Palmé A.E., Cheddadi R., Latta R.G. 2004 . Impact of Ice Ages on the genetic structure of trees and shrubs. Philos. Trans. R. Soc. Lond. B Biol. Sci. 359 : 197 – 207 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Latta R.G., Linhart Y.B., Fleck D., Elliot M. 1998 . Direct and indirect estimates of seed versus pollen movement within a population of ponderosa pine. Evolution 52 : 61 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Latta R.G., Mitton J.B. 1999 . Historical separation and present gene flow through a zone of secondary contact in ponderosa pine. Evolution 53 : 769 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lesser M.R., Jackson S.T. 2012 . Making a stand: five centuries of population growth in colonizing populations of Pinus ponderosa. Ecology 93 : 1071 – 1081 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Lesser M.R., Jackson S.T. 2013 . Contributions of long-distance dispersal to population growth in colonising Pinus ponderosa populations. Ecol. Lett. 16 : 380 – 389 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Liaw A., Wiener M. 2002 . Classification and regression by randomForest. R News 2 : 18 – 22 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Little E.L. 1971 . Atlas of United States Trees, vol. 1 . Conifers and important hardwoods. Washington (DC) : USDA Forest Service. Miscellaneous Publication 1146. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Maguire K.C., Nieto-Lugilde D., Blois J.L., Fitzpatrick M.C., Williams J.W., Ferrier S., Lorenz D.J. 2016 . Controlled comparison of species- and community-level models across novel climates and communities. Proc. R. Soc. Lond. B Biol. Sci. 283 : 20152817 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Marcer A., Méndez-Vigo B., Alonso-Blanco C., Picó F.X. 2016 . Tackling intraspecific genetic structure in distribution models better reflects species geographical range. Ecol. Evol. 6 : 2084 – 2097 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC McCune B., and Mefford J. 2009 . HyperNiche. Nonparametric multiplicative habitat modeling, Version 2.28. Gleneden Beach (OR) : MjM Software. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Milborrow S. 2016 . Earth: multivariate adaptive regression splines. R package version 4.4.4. Available from: URL https://CRAN.R-project.org/package=earth. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Namkoong G., Conkle M.T. 1976 . Time trends in genetic control of height growth in ponderosa pine. Forest Sci. 22 : 2 – 12 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Neale D.B., Sederoff R.R. 1989 . Paternal inheritance of chloroplast DNA and maternal inheritance of mitochondrial DAN in loblolly pine. Theor. Appl. Genet. 77 : 212 – 216 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Niebling C.R., Conkle M.T. 1990 . Diversity of Washoe pine and comparisons with allozymes of ponderosa pine races. Can. J. Forest Res. 20 : 298 – 308 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Norris J.R., Jackson S.T., Betancourt J.L. 2006 . Classification tree and minimum-volume ellipsoid analyses of the distribution of ponderosa pine in the western USA. J. Biogeogr. 33 : 342 – 360 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Norris J.R., Betancourt J.L., Jackson S.T. 2016 . Late Holocene expansion of ponderosa pine (Pinus ponderosa) in the Central Rocky Mountains, USA. J. Biogeogr. 43 : 778 – 790 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Oliver W.W., Ryker R.A. 1990 . Ponderosa pine. In: Burns R.M., Honkala B.H., editors. Silvics of North America: 1. Conifers, vol 1 . Agricultural Handbook 654. Washington : U.S. Department of Agriculture Forest Service. p. 413 – 424 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Oney B., Reineking B., O’Neill G., Kreyling J. 2013 . Intraspecific variation buffers projected climate change impacts on Pinus contorta. Ecol. Evol. 3 : 437 – 449 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pauls S.U., Nowak C., Bálint M., Pfenninger M. 2013 . The impact of global climate change on genetic diversity within populations and species. Mol. Ecol. 22 : 925 – 946 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pearman P.B., D’Amen M., Graham C.H., Thuiller W., Zimmermann N.E. 2010 . Within-taxon niche structure: niche conservatism, divergence and predicted effects of climate change. Ecography 33 : 990 – 1003 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Petit R.J., Kremer A., Wagner B. 1993 . Finite island model for organelle and nuclear genes in plants. Heredity 71 : 630 – 641 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pfenninger M., Nowak C., Magnin F. 2007 . Intraspecific range dynamics and niche evolution in Candidula land snail species. Biol. J. Linn. Soc. 90 : 303 – 317 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Potter K.M., Hipkins V.D., Mahalovich M.F., Means R.E. 2013 . Mitochondrial DNA haplotype distribution patterns in Pinus ponderosa (Pinaceae): range-wide evolutionary history and implications for conservation. Am. J. Bot. 100 : 1562 – 1579 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Potter K.M., Hipkins V.D., Mahalovich M.F., Means R.E. 2015 . Nuclear genetic variation across the range of ponderosa pine (Pinus ponderosa): phylogeographic, taxonomic and conservation implications. Tree Genet. Genomes 11 : 38 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Potter K.M., Jetton R.M., Bower A., Jacobs D.F., Man G., Hipkins V.D., Westwood M. 2017 . Banking on the future: progress, challenges and opportunities for the genetic conservation of forest trees. New Forests 48 : 153 – 180 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Prasad A.M., Potter K.M. 2017 . Macro-scale assessment of demographic and environmental variation within genetically derived evolutionary lineages of eastern hemlock (Tsuga canadensis), an imperiled conifer of the eastern United States. Biodivers. Conserv. 10.1007/s10531-017-1354-4 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref R Core Team 2014 . R: a language and environment for statistical computing. Vienna (Austria) : R Foundation for Statistical Computing. Available from: URL http://www.R-project.org/. Rehfeldt G.E. 1986b . Adaptive variation in Pinus ponderosa from Intermountain regions. II. Middle Columbia River System. Research Paper INT-RP-373. Ogden, UT : U.S. Department of Agriculture, Forest Service, Intermountain Research Station. p9 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E. 1986b . Adaptive variation in Pinus ponderosa from intermountain regions. II. Middle Columbia river system. Research Paper INT-RP-373. Ogden, UT : U.S. Department of Agriculture, Forest Service, Intermountain Research Station. p9 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E. 1990 . Genetic differentiation among populations of Pinus ponderosa from the Upper Colorado River Basin. Bot. Gaz. 151 : 125 – 137 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E. 2006 . A spline model of climate for the Western United States. General Technical Report RMRS-GTR-165. Fort Collins, CO : U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station. p.21 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rehfeldt G.E., Jaquish B.C., López-Upton J., Sáenz-Romero C., St Clair J.B., Leites L.P., Joyce D.G. 2014 . Comparative genetic responses to climate for the varieties of Pinus ponderosa and Pseudotsuga menziesii: realized climate niches. Forest Ecol. Manag. 324 : 126 – 137 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rosauer D.F., Catullo R.A., VanDerWal J., Moussalli A., Moritz C. 2015 . Lineage range estimation method reveals fine-scale endemism linked to Pleistocene stability in Australian rainforest herpetofauna. PLoS One 10 : e0126274 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Roberts D.R., Hamann A. 2016 . Climate refugia and migration requirements in complex landscapes. Ecography 39 : 1238 – 1246 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Rudloff E., von Lapp M.S. 1992 . Chemosystematic studies in the genus Pinus. VII. The leaf oil terpene composition of ponderosa pine, Pinus ponderosa. Can. J. Bot. 70 : 374 – 378 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Shinneman D.J., Means R.E., Potter K.M., Hipkins V.D. 2016 . Exploring climate niches of ponderosa pine (Pinus ponderosa Douglas ex Lawson) haplotypes in the western United States: implications for evolutionary history and conservation. PLoS One 11 : e0151811 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Schoener T.W. 1970 . Nonsynchronous spatial overlap of lizards in patchy habitats. Ecology 51 : 408 – 418 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sorensen F.C., Mandel N.L., Aagaard J.E. 2001 . Role of selection versus historical isolation in racial differentiation of ponderosa pine in southern Oregon: an investigation of alternative hypotheses. Can. J. Forest Res. 31 : 1127 – 1139 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sorensen F.C., Weber J.C. 1994 . Genetic variation and seed transfer guidelines for ponderosa pine in the Ochoco and Malheur National Forests of central Oregon. USDA Forest Service Research Paper PNW-RP-472, Washington DC , USA . 24 p. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Sork V.L., Davis F.W., Westfall R., Flint A., Ikegami M., Wang H., Grivet D. 2010 . Gene movement and genetic association with regional climate gradients in California valley oak (Quercus lobata Née) in the face of climate change. Mol. Ecol. 19 : 3806 – 3823 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Valladares F., Matesanz S., Guilhaumon F., Araújo M.B., Balaguer L., Benito-Garzón M., Cornwell W., Gianoli E., van Kleunen M., Naya D.E., Nicotra A.B., Poorter H., Zavala M.A. 2014 . The effects of phenotypic plasticity and local adaptation on forecasts of species range shifts under climate change. Ecol. Lett. 17 : 1351 – 1364 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wang C.-W. 1977 . Genetics of ponderosa pine. Washington (DC) : U.S. Department of Agriculture, Forest Service . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wells O.O. 1964 . Geographic variation in Ponderosa Pine. Silvae Genetica 13 : 89 – 103 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Williams S.E., Shoo L.P., Isaac J.L., Hoffmann A.A., Langham G. 2008 . Towards an integrated framework for assessing the vulnerability of species to climate change. PLoS Biol. 6 : e325 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wisz M.S., Pottier J., Kissling W.D., Pellissier L., Lenoir J., Damgaard C.F., Dormann C.F., Forchhammer M.C., Grytnes J.-A., Guisan A., Heikkinen R.K., Høye T.T., Luoto M., Maiorano L., Nilssson M.-C., Normand S., Öckinger E., Schmidt N.M., Termansen M., Timmerman A., Wardle D.A., Aastrup P., Svenning J.-C. 2013 . The role of biotic interactions in shaping distributions and realized assemblages of species: implications for species distribution modelling. Biol. Rev. 88 : 15 – 30 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wofford A.M., Finch K., Bigott A., Willyard A. 2014 . A set of plastid loci for use in multiplex fragment length genotyping for intraspecific variation in Pinus (Pinaceae). Appl. Plant Sci. 2 : 1400002 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Yannic G., Pellissier L., Ortego J., Lecomte N., Couturier S., Cuyler C., Dussault C., Hundertmark K.J., Irvine R.J., Jenkins D.A., Kolpashikov L., Mager K., Musiani M., Parker K.L., Røed K.H., Sipko T., Þórisson S.G., Weckworth B.V., Guisan A., Bernatchez L., Côté S.D. 2014 . Genetic diversity in caribou linked to past and future climate change. Nat. Clim. Change 4 : 132 – 137 . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Published by Oxford University Press, on behalf of the Society of Systematic Biologists 2018. This work is written by (a) US Government employee(s) and is in the public domain in the US. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) Published by Oxford University Press, on behalf of the Society of Systematic Biologists 2018. This work is written by (a) US Government employee(s) and is in the public domain in the US.

Journal

Systematic Biology – Oxford University Press

Published: Nov 1, 2018

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Intraspecific Niche Models for Ponderosa Pine (Pinus ponderosa) Suggest Potential Variability in Population-Level Response to Climate Change

Intraspecific Niche Models for Ponderosa Pine (Pinus ponderosa) Suggest Potential Variability in Population-Level Response to Climate Change

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Intraspecific Niche Models for Ponderosa Pine (Pinus ponderosa) Suggest Potential Variability in Population-Level Response to Climate Change

Intraspecific Niche Models for Ponderosa Pine (Pinus ponderosa) Suggest Potential Variability in Population-Level Response to Climate Change

References (121)

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies