Intraspecific variation buffers projected climate change impacts on Pinus contorta.

1.1 MB·English

by  Oney, Brian

Intraspecific variation buffers projected climate change
Pinus contorta
impacts on
Brian Oney1, Björn Reineking1, Gregory O’Neill2 & Juergen Kreyling3
1BiogeographicalModelling,UniversityofBayreuth,BayCEER,Bayreuth,Germany
2Research,InnovationandKnowledgementManagementBranch,BCMinistryofForestsandRange,Vernon,BritishColumbia,Canada
3Biogeography,UniversityofBayreuth,BayCEER,Bayreuth,Germany
Keywords Abstract
Intraspecificdiversity,lodgepolepine,niche
Species distribution modeling (SDM) is an important tool to assess the impact
modeling,NorthAmerica,rangeshift,
within-speciesvariability. of global environmental change. Many species exhibit ecologically relevant
intraspecific variation, and few studies have analyzed its relevance for SDM.
Correspondence Here, we compared three SDM techniques for the highly variable species Pinus
JuergenKreyling,Biogeography contorta. First, applying a conventional SDM approach, we used MaxEnt to
UniversityofBayreuth,BayCEER model the subject as a single species (species model), based on presence–
Universita¨tsstr,40Bayreuth,
absence observations. Second, we used MaxEnt to model each of the three most
Germany95447.Tel:[phone];
prevalent subspecies independently and combined their projected distributions
Fax:[phone];
E-mail:[email protected] (subspecies model). Finally, we used a universal growth transfer function
(UTF), an approach to incorporate intraspecific variation utilizing provenance
FundingInformation trial tree growth data. Different model approaches performed similarly when
ThispublicationwasfundedbytheGerman predicting current distributions. MaxEnt model discrimination was greater
ResearchFoundation(DFG)andthe (AUC – species model: 0.94, subspecies model: 0.95, UTF: 0.89), but the UTF
UniversityofBayreuthinthefunding was better calibrated (slope and bias – species model: 1.31 and (cid:1)0.58, subspe-
programmeOpenAccessPublishing.
cies model: 1.44 and (cid:1)0.43, UTF: 1.01 and 0.04, respectively). Contrastingly,
Received:11October2012; Accepted:15 for future climatic conditions, projections of lodgepole pine habitat suitability
October2012 diverged. In particular, when the species’ intraspecific variability was acknow-
ledged, the species was projected to better tolerate climatic change as related to
EcologyandEvolution2013;3(2):437–449 suitable habitat without migration (subspecies model: 26% habitat loss or UTF:
24% habitat loss vs. species model: 60% habitat loss), and given unlimited
doi:10.1002/ece3.426
migration may increase amount of suitable habitat (subspecies model: 8% habi-
tat gain or UTF: 12% habitat gain vs. species model: 51% habitat loss) in the
climatic period 2070–2100 (SRES A2 scenario, HADCM3). We conclude that
models derived from within-species data produce different and better projec-
tions, and coincide with ecological theory. Furthermore, we conclude that
intraspecific variation may buffer against adverse effects of climate change.
A key future research challenge lies in assessing the extent to which species can
utilize intraspecific variation under rapid environmental change.
tedtoextinctionintheclimaticperiod2070–2100ifunable
Introduction
todisperse(Thuilleret al.2005).Speciesdistributionmod-
Projections of climate change and the related impacts on eling (SDM) encompasses a broad range of techniques
species distributions suggest significant ecological distur- (GuisanandZimmermann2000),andapplicationsranging
bance, especially when considering species’ range losses from paleobiology (Svenning et al. 2011) to spread and
(e.g.,Thomaset al.2004;Urbanet al.2012).Modelingthe control of infectious vector-borne disease (Fischer et al.
future potential distribution of 1,350 European plant 2011). However, major uncertainties are associated with
species under various greenhouse-gas emission scenarios SDM, such as the realistic modeling of migration rates
indicatesthatmorethanhalfofthesespeciescouldbecome (Best et al. 2007; Nathan et al. 2011), biotic interactions
vulnerable, endangered, critically endangered, or commit- (Araújo and Luoto 2007; Preston et al. 2008; Meier et al.
ª2013TheAuthors.PublishedbyBlackwellPublishingLtd.ThisisanopenaccessarticleunderthetermsoftheCreative 437
CommonsAttributionLicense,whichpermitsuse,distributionandreproductioninanymedium,provided
theoriginalworkisproperlycited.
ImportanceofIntraspecificVariationinSDM B.Oneyetal.
2010),forexample,fecundityandcompetition(Clarket al. importance (Krieger 2001; Karst 2010). Naturally and
2011), consideration of micro-climate (Dobrowski 2011; through human influence, the economically important
Hof et al. 2011; Suggitt et al. 2011) or climate extremes subspecieslatifoliahasbeenadvancingnorthwardsincethe
(Zimmermann et al. 2009), and persistence of ecosystem end of the last glacial maximum (MacDonald 1991; Faze-
structure and reliance upon other predictions (Dormann kas and Yeh 2006) and recent observations hint toward a
2007;Wienset al.2009).Nonetheless,successfulretrospec- continued northward spread (Johnstone and Chapin
tive predictions of shifts in bird population sizes demon- 2003). Pinus contorta exhibits great ecological variation,
stratethevalueofspeciesdistributionmodels(Greenet al. andcanbeseparatedintothreegeneticallyandecologically
2008). Currently, migration rates and pathways are added distinct subspecies: contorta, murrayana and latifolia
to the models, thereby increasing realism of the projected (LotanandCritchfield1990),makingitanidealcandidateto
resultswhilestillbeingrestrictedbylimitedspecies-specific investigatetheeffectsofintraspecificvariabilityonSDM.
knowledge on potential migration rates (Fischer et al. Wearguethattheintegrationofwithin-speciesvariation
2011). Austin (2007) presented a general framework for is a necessary step when considering ecological theory in
SDMstudies,suchthatecologicaltheorymustconcurwith SDM.Weinvestigatedtheeffectofincorporatingintraspe-
a data model and statistical model. The incorporation of cific variation on SDM performance and projections using
variationwithinaspeciesorcloselyrelatedgroupofspecies range-wide, geo-referenced, subspecies information and
has been recently addressed as an important challenge to provenance test data of Pinus contorta. Two approaches
improve SDM (Zimmermann et al. 2010). However, “spe- driven by intraspecific data were compared with conven-
cies” is a taxonomic designation, and may not necessarily tional SDM, which allowed for robust comparison of our
designate an ecologically homogeneous group of organ- findings against a particular way to incorporate intraspe-
isms, especially when intraspecific ecotypes occur. Experi- cific variability. We expected a difference in modeling
mental evidence suggests that conventional SDM cannot results between approaches that ignore intraspecific varia-
capturetheclimaticresponseofspeciesbytreatingthemas tion and those that incorporate intraspecific variation.
homogeneous units (Beierkuhnlein et al. 2011). However, Furthermore, given the different theoretical backgrounds
O’Neill et al. (2008) present an interesting approach to of the modeling approaches, we expected substantial
account for intraspecific variation of growth using prove- differencesinperformancebetweenmodelingapproaches.
nance trial data. Similarly, Benito Garzo´n et al. (2011)
investigated future tree survival when considering popula-
Methods
tionvariation,andfoundthatacknowledgmentofthisvari-
ability provides a more positive projected outlook into
Climate data
future climates when compared with conventional SDM
approaches. Furthermore, Pearman et al. (2010) note a The worldclim climate data (Hijmans et al. 2005) was
slight improvement of their models when incorporating used with a resolution of 2.5 arcminutes, which was
within-species variation when compared with the tradi- reprojected to Albers Equal Area projection (resolution of
tionalSDMmethod.HamannandWang(2006)developed 4 km2) with GRASS GIS (GRASS Development Team
an interesting community modeling approach to account 2012). Climate variables were derived using the worldclim
for ecosystem level variability, which has been successfully dataset. Climate variable derivation formulae (Table 1)
applied to highly prevalent tree populations, including were taken from two sources: those used in O’Neill et al.
Pinus contorta, in order to consider future development of (2008), and additional bioclimatic variables from Hijmans
local forest stands under climate change scenarios (Gray et al. (2005). Present climate data (1950–2000) were used
and Hamann 2012) as well as past climates (Roberts and for model training and validation, and future climate sim-
Hamann 2012). The “ecosystem-based” approach has sev- ulations were used for SDM model projections (HadCM3,
eraladvantagesoverspeciesmodels,onebeingtheabilityto A2a emissions scenario; Nakićenović et al. 2000) for peri-
model individual populations in a changing environment ods 2010–2040, 2040–2070, and 2070–2100. Considering
(Roberts and Hamann 2012). Within-species variation, current CO2 emissions as well as plausible future political
according to the insurance hypothesis (Yachi and Loreau and social developments, the use of A2 simulation data
1999), contributes to a species’ ability to utilize various appears justified (Moss et al. 2010).
resources (Joshi et al. 2001; Kreyling et al. 2012a,b) and
therebyadapttoecologicalchange(Daviset al.2005;Nus-
Occurrence data
seyet al.2005;Skellyet al.2007).
Pinus contorta is a pyrophilic, widely distributed, Occurrence data (Fig. 1) were obtained mainly from the
outcrossing, wind-dispersed conifer species of high Vegetation Resource Inventory of the British Columbia
ecological (Lotan and Critchfield 1990) and economic (BC) Forest Service, provided at a 1600-m grid and the
438 ª2013TheAuthors.PublishedbyBlackwellPublishingLtd.
B.Oneyetal. ImportanceofIntraspecificVariationinSDM
Table1. Climate variables investigated across the ranges of Pinus
contortaandsubspecies. Modeling
Range
Climatevariables Unit Source
AnnualHeat/MoistureIndex ◦C/mm Wangetal.(2006)
SummerHeat/MoistureIndex ◦C/mm Wangetal.(2006)
MeanAnnualTemp. ◦C Wangetal.(2006)
MeanWarmMonthlyTemp. ◦C Wangetal.(2006)
MeanColdMonthlyTemp. ◦C Wangetal.(2006)
Temp.Difference/AnnualRange ◦C Wangetal.(2006)
MeanAnnualPrecip. mm Wangetal.(2006)
MeanSummerPrecip. mm Wangetal.(2006)
Isothermality ◦C/◦C Hijmansetal.(2005)
MeanDiurnalRange ◦C Hijmansetal.(2005)
Temp.Seasonality ◦C Hijmansetal.(2005)
MeanTemp.ofWettestQuarter ◦C Hijmansetal.(2005)
MeanTemp.ofDriestQuarter ◦C Hijmansetal.(2005)
MeanTemp.ofWarmestQuarter ◦C Hijmansetal.(2005)
MeanTemp.ofColdestQuarter ◦C Hijmansetal.(2005) contorta
Precip.ofWettestMonth mm Hijmansetal.(2005) murrayana
Precip.ofDriestMonth mm Hijmansetal.(2005) latifolia
Precip.Seasonality mm Hijmansetal.(2005)
Precip.ofWettestQuarter mm Hijmansetal.(2005)
Figure1. Pinus contorta subspecies distributions across its natural
Precip.ofDriestQuarter mm Hijmansetal.(2005)
range. The number of observations documented to subspecies was
Precip.ofColdestQuarter mm Hijmansetal.(2005)
low; therefore the nearest observation that was classified to
Precip.ofWarmestQuarter mm Hijmansetal.(2005)
MeanMax.Temp.ofDriestQuarter ◦C Thisstudy subspecies was assigned that observations’ subspecies classification.
The dashed line outlines the model building buffer area i.e. model
trainingandevaluation(seeSupportingInformation).Occurrencedata
were obtained mainly from the British Columbia Ministry of Forests
US Forest Service, provided on a 10-km grid. Also, occur- andUSForestServiceForestInventoryandAnalysis;theobservational
rence data were further supplemented with online data cover most of the natural range. Within the modeling range,
resources such as herbaria, botanical gardens, and plant eachrastercellwasassignedeitherapresenceorabsence.
databases (see Supporting Information). Each dataset was
examined extensively for outliers, relative to the most observations does not equal the total observations due to
comprehensive P. contorta occurrence data from Little co-occurrence (see Supporting Information).
(1971). Unfortunately, most occurrences were not docu-
mented and classified to subspecies. In order to address
Range designation
this problem, occurrences without a subspecies classifica-
tion were assigned the subspecies of their nearest spatial In order to avoid biased niche estimates arising from
neighbor (see Supporting Information). The nearest treating areas as absence that are climatically suitable but
neighbor analysis assigned the subspecies contorta to some are not occupied because of geographic isolation, a buffer
observations that had elevations far outside its observed of 1,000 km around the occurrence dataset was calculated
elevational range. Lotan and Critchfield (1990) describe for the entire species as well as for each subspecies for
P. contorta subspecies contorta as occurring mainly model building and evaluation. Each model was trained
between sea-level and 610-m altitude. One of the P. con- on its own buffer (see Supporting Information), evaluated
torta subspecies contorta provenances used in the on probabilities within the whole species buffer (Fig. 1 –
Illingworth trial (O’Neill et al. 2008) was sampled at an “Modeling Range”), and projected to the North American
elevation of 1266 m (sampled with the reprojected continent north of 23°N. Pinus contorta has yet to be doc-
worldclim-SRTM elevation database); we used this value umented to occur north of 65°N or in Saskatchewan,
as the upper elevational limit of P. contorta subspecies except in the southwest corner (Little 1971). We were
contorta. These observations were in areas in the Pacific unable to obtain a thorough dataset for Alberta, whose
Northwest where subspecies latifolia and subspecies area was, except for presence data, omitted from the anal-
contorta introgress, and were therefore assigned to subspe- yses. Areas lying within the buffer of 1000 km, but well
cies latifolia. Observations within a raster cell for each outside the documented natural range of P. contorta were
subspecies contorta: 2048, murrayana: 1449, latifolia: included as absence data. This meant that, within the buf-
42,342, and totaled: 45,785. The sum of the subspecies fer, areas north of 65°N, west of the Alberta-Saskatchewan
ª2013TheAuthors.PublishedbyBlackwellPublishingLtd. 439
ImportanceofIntraspecificVariationinSDM B.Oneyetal.
border, or areas much further (>500 km) away from the model probabilities P at a given location x were then
subsp,i
edges of the observations were considered to be absence combined (Eq. 1) as in Pearman et al. (2010) to yield the
data. Areas of southern Alaska were omitted from analysis occurrence probability of the species as a whole, P :
tot
onthegroundsthatthesamplingintensityissparseinthis YN
region, and because two observations (Arctos Museum of PtotðxÞ¼1(cid:1) ð1(cid:1)Psubsp;iðxÞÞ (1)
the North, 2010) indicate possibly viable populations in i¼1
south-centralAlaska(seeSupportingInformation).
Universal transfer function approach
Modeling approaches
The Illingworth provenance trial (Illingworth 1978) began
ThreeSDMapproacheswereusedtomodeldistributionsof in the late 1960′s, taking seed from 140 provenances
P. contorta. MaxEnt(Phillips et al.2006) was used (1)as a across most of the range of P. contorta, growing seeds in
conventionalSDM(GuisanandZimmermann2000;Wiens
a nursery, and planting the 3-year-old seedlings across BC
et al.2009),thespeciesobservationaldataweremodeledas
and Yukon Territory in 1974. An incomplete factorial
if the species was ecotypically homogeneous, called here-
design was used to test 60 populations at each of 60 test
after the “species model”. In order to incorporate intra-
sites. We derived the UTF according to the same methods
specificvariation,themostprevalentP. contortasubspecies
from O’Neill et al. (2008), using the Illingworth trial
(contorta, murrayana, and latifolia) were modeled; (2) as
35- year-old tree plot growth data (volume per hectare [V
autonomous units, call hereafter the “subspecies model”. PH], m3/ha), a measure that combines height, diameter
Finally, the “Universal Transfer Function” (UTF) from
and survival. Using the original model structure, tree plot
O’Neill et al. (2008) was used (3) to incorporate the
growth data were used to fit the UTF to the worldclim
observedvariabilityinprovenancetrialtreegrowth.
climate variables. The UTF projections were produced in
two steps. First, at each test, site (S), individual prove-
Species model nance (P), mean cold month temperature (MCMTP) were
fitted to population production (VPH ) using a unimodal
P
Candidate climate variables were tested for collinearity
Cauchy function to develop an individual transfer func-
with each other with Spearman’s non- parametric correla-
tion for each test site:
tion. Correlation among candidate climate variables was A
examined. Where pairs of variables were highly correlated VPHP ¼ (cid:2) S (cid:3) (2)
(q > 0.7), a univariate generalized additive model (GAM) 1þ MCMCTP(cid:1)BS 2
S
was fitted to the test data (see below) using each highly
where A , B and C are the fitted function parameters,
correlated variable. In order to obtain less correlated vari- S S S
yielding a transfer function per trial site S (N = 42, see
ables, the variable of each pair that yielded the greater
Eq. 2; for background see Raymond and Lindgren (1990);
AIC was omitted. MaxEnt (Phillips et al. 2006; Maximum
Lindgren and Ying (2000)). Inter-site variation in the
Entropy, version 3.3a) was then used to model the cur-
transfer functions, that is, the variation in the transfer
rent distribution with the full occurrence dataset. Austin
function parameters A , B and C was then modeled as a
(2007) discussed species response curves that coincide S S S
function of site climate, using mean annual temperature,
with ecological theory and we assumed a smooth response
mean coldest month temperature, and annual heat/mois-
along a climatic gradient. Therefore, MaxEnt “feature”
ture index (see O’Neill et al. 2008; for details). The result-
types linear, product and quadratic features were used
ing UTF (Eq. 2) predicts 35-year-old tree growth of a
(Elith et al. 2011), and the “samples-with-data” (SWD)
population from any climate growing in any climate
input data format (see Appendix S4 in Elith et al. 2011).
(VPH ), and was used to predict current and future
Otherwise, default settings were used. We use the term S
growth estimates. The growth estimates for each raster
“absence data” to refer to the MaxEnt analogy “back-
cell were then used to explain the observational data
ground samples”, and model absence data as background
using logistic regression. The logistic regression model
psamples (MaxEnt argument: environmentallayers).
was then used to convert present and future growth esti-
mates to occurrence probabilities.
Subspecies model
Each subspecies was treated as its own viable species, and
Model evaluation
was individually investigated for determinant climate vari-
ables. This was done using the same algorithm as for the Observations were randomly partitioned into two data-
“species model” approach. The resultant N subspecies sets, with 70% of the cells selected for training and the
440 ª2013TheAuthors.PublishedbyBlackwellPublishingLtd.
B.Oneyetal. ImportanceofIntraspecificVariationinSDM
remaining 30% of the cells used for model evaluation, as we did not set a probability threshold to create a presence
is common in SDM studies (e.g., Thuiller et al. 2004; –absence map from the model projections, we took the
Fischer et al. 2011). The same set of cells was used for occurrence data and quantitated the area of suitable
training and evaluating all models, that is, across all three habitat at each occurrence location. All data preparation,
approaches. Model performance was assessed by the area analyses, and visualization were conducted in R (R Devel-
under the receiving operator characteristic curve (AUC), opment Core Team 2012) and GRASS GIS (GRASS
which measures discrimination, the Nagelkerke-Cox- Development Team 2012, version 6.4).
Snell-Maddala-Magee R2–index, a general indicator of
performance, and the logistic calibration curve, of which
Results
the slope indicates over- or underfitting (calibration) and
the intercept indicates bias (see Reineking and Schro¨der
Model performance
2006). A calibration curve slope greater than 1 (a > 1)
indicates under-fitting/over-regularization, vice versa for For each MaxEnt model target, a separate set of climate
a < 1, and bias indicates how well the numerical values variables was found to be best suited for describing the
of the actual and predicted probabilities correspond; that occurrence distributions, and the response curve charac-
is, an intercept of less than 1 (b > 1) indicates that the teristics give insight into the relationship between
predicted probabilities are too high and vice versa. The P. contorta distribution and climate (Table 2). The shapes
variable “importance” was calculated by randomizing of the curves representing the response of each subspecies
(n = 10) each variable and then subtracting the mean provide more evidence that each subspecies should react
change in AUC from the full model, and each AUC differently to climatic changes. It is interesting to note the
difference is divided by the sum of all AUC changes to variety with which the models depict each subspecies,
derive the relative importance. especially when comparing with the “species” model. In
general, the models performed well. According to the R2
-index and AUC values, the species and subspecies models
Suitable area quantitation performed generally better than the UTF and were more
discriminative, respectively (Fig. 2). According to the cali-
In order to avoid the somewhat arbitrary procedure of
bration analyses, the occurrence probabilities of the spe-
threshold choice (Fielding and Bell 1997), the total area
cies and subspecies models were under-fitted (slope: 1.31
of habitat suitability Atot was quantitated by multiplying and 1.44) and overestimated (bias: (cid:1)0.58 and (cid:1)0.43),
the predicted probability of occurrence Pj for each cell j
respectively, while the predicted occurrence probabilities
with the area of the cell Areaj and summing over the M
from the UTF are well calibrated (slope: 1.01 and bias:
analyzed cells (Eq. 3).
1.44).
XM
It is notable that the species model predicts low occur-
Areatot ¼ Pj(cid:3)Areaj (3) rence probabilities of subspecies murrayana in the Sierra
j¼1
Nevada. All models predict a medium to high occurrence
The AUC is a rank-based measure, which allowed us to probability in eastern Canada and southwest Alaska where
quantitate the discriminative ability of the models, with- it is absent. The MaxEnt species model predicts a high
out choosing a threshold and creating a projected binary occurrence probability in the Northwest Territories, where
map of presences and absences. We refer to habitat suit- P. contorta is currently extending its range, but not to the
ability as the modeled probability of occurrence, which is extent predicted.
important for predicting the success of seedlings, but also
serves as reference for production potential, although
Importance of intraspecific variability for
adult trees have been observed to tolerate a wide range of
future projections
conditions. The relative amounts of suitable habitat were
quantitated with two dispersal scenarios: “no dispersal” Predicted future habitat suitability shows a general north-
and “full dispersal”. A full dispersal scenario is the ward shift of suitable habitat (compare Figs. 2 & 3). The
optimistic extreme and assumes that the species can reach models predict different future development of suitable
all habitats on the North American continent. The no area, with a marked difference between the UTF or the
dispersal scenario is the pessimistic extreme, and assumes subspecies model and the species model (Fig. 4); that is,
that the species will not be able to migrate. These dis- the difference is between the models that do or do not
persal scenarios were used to bound and illustrate the incorporate intraspecific variability. More specifically, in a
range of possible future locations of the species’ climate no-dispersal scenario, the UTF predicts 24% habitat loss
habitat, and to guide assisted migration efforts. Because for the period of 2070–2100 relative to the current
ª2013TheAuthors.PublishedbyBlackwellPublishingLtd. 441
ImportanceofIntraspecificVariationinSDM B.Oneyetal.
Table2. ThesummarizingdimensionsoftheclimatenicheforeachmodelingsubjectusingMaxEnt.
Subject Variable I.(%) Range Curve
Subsp.contorta Meandiurnalrange 44.7 4–13 Unimodal,pos.
Maximumtemperatureofdryquarter 29.8 (cid:1)7–30 Unimodal,pos.
Meantemperatureofwetquarter 25.5 (cid:1)7–19 Unimodal,pos.
Subsp.murrayana Precipitationofwarmquarter 40.3 1–232 Sigmoidal,neg.
Precipitationofcoldquarter 25.3 11–581 Sigmoidal,neg.
Maximumtemperatureofdryquarter 14.9 (cid:1)7–21 Unimodal,pos.
Precipitationseasonality 9.8 19–110 Unimodal,pos.
Precipitationofwetquarter 9.7 152–760 Sigmoidal,neg.
Subsp.latifolia Meanannualtemperature 72.8 (cid:1)12–14 Unimodal,pos.
Precipitationofdryquarter 14.7 0–262 Unimodal,pos.
Precipitationseasonality 12.5 6–70 Sigmoidal,neg.
Speciesmodel Meantemperatureofwarmquarter 72.8 0–22 Unimodal,pos.
Precipitationofdryquarter 19.8 2–470 Unimodal,pos.
Precipitationseasonality 7.4 20–75 Sigmoidal,neg.
Asetofvariableswasfoundtobelesscorrelated(q<0.7)andimportanttoexplainthedistributionsofeachsubspecies.Therangesofvariables
(unitsareinTable1),inwhicheachsubspeciesrespondstoclimaticconditions,andtheshapeoftheirresponsecurvehighlightthedifferentchar-
acteristics of each subspecies, and provide an informative contrast to the conventional method (species model). Furthermore, the importance (I.
(%)) of each variable, measured by the relative change in AUC when singularly randomizing, hints at the dominant climate factors depicted by
eachmodel.
scenario (1950–2000), and similarly the subspecies model than either the species or subspecies models, whereas the
predicts 26% habitat loss, while the species model species and subspecies models were more discriminative,
predicts 60% habitat loss. Furthermore, given perfect but were poorly calibrated. The UTF approach used the
dispersal to newly suitable climatic areas for the period of occurrence data only via logistic regression to convert the
2070–2100, the UTF predicts 12% habitat gain, and the predicted tree growth to occurrence probabilities;
subspecies model similarly predicts 8% habitat gain, but independently, the growth model was derived from the
the species model predicts 51% habitat loss. As can be Illingworth provenance trial data. While the UTF
seen in Fig. 4, the predicted change in climatically suit- approach still discriminates well, we see a non-trivial per-
able area is gradual over time. Examining the subspecies formance difference between the species or subspecies
model projections more carefully (Fig. 5), predicted models, and the UTF. Notably, the subspecies model
climatic change affects the subspecies differently: in a no- slightly outperformed the species model, coinciding with
dispersal scenario, subspecies latifolia may lose current the results from Pearman et al. (2010).
habitat (23% habitat loss) by the period 2070–2100, and The SDMmodelalgorithm (e.g.,RandomForest,GAM,
subspecies contorta is also predicted to experience habitat MaxEnt, etc.) has been shown to be an important source
loss (28% habitat loss), whereas the subspecies murrayana of uncertainty (Dormann et al. 2008). However, in this
is predicted to drastically lose current habitat suitability study, we observed that incorporating intraspecific varia-
(95% habitat loss). In a full dispersal scenario, subspecies tion is another important source of variation in addition
latifolia may experience habitat gain (12% habitat gain), to model algorithm. Results of the species model and the
whereas subspecies contorta is predicted to experience a other two models differed greatly in terms of predicted
decrease in habitat suitability (95% habitat loss) as well as changes in the extent of future suitable habitat. In con-
subspecies murrayana (95% habitat loss). trast, we did not observe a notable difference between the
approaches driven by intraspecific data, that is, the
subspecies model and the UTF (Fig. 4). In the two
Discussion
dispersal scenarios, the intraspecific approaches agreed
The three models exhibited similar performance and with each other, but differed from the conventional spe-
performed well for current distributions and climates. cies model, which predicted much greater climate change
The species model was more discriminative than the effects. Given that the two intraspecific approaches come
UTF, but was under-fitted and yielded overestimated to similar conclusions concerning climate change impacts
occurrence probabilities similar to the subspecies model. suggests that the result does not depend on the particular
The UTF was less discriminative but better calibrated way in which intraspecific variation is taken into account,
442 ª2013TheAuthors.PublishedbyBlackwellPublishingLtd.
B.Oneyetal. ImportanceofIntraspecificVariationinSDM
(a) (a)
(b) (b)
(c) (c)
Figure2. Occurrence probability distributions were modeled in the Figure3. Projections of the three models (a) species model (b)
current climatic period (1950–2000) and the three models (a) species subspecies model and (c) universal transfer function onto North
model (b) subspecies model and (c) universal transfer function (UTF) Americanclimatefortheperiod2070–2100undertheA2aemissions
were projected onto North American climate data (Hijmans et al. scenariofromtheHadCM3globalclimatemodel.
2005). The models were trained and evaluated on the modeling
range(seeFig.1)andprojectedtotheprojectionrange(dashedline), P. contorta climatically suitable habitat within BC from
which corresponds to most of the North American continent. The Hamann and Wang (2006) as well as Gray and Hamann
MaxEnt models are highly discriminative, but the UTF is better
(2012), whose community modeling approach yielded
calibrated.
promising results for projecting ecolimatic shifts in forest
and thereby corroborates previous research using the stands (Gray et al. 2011; Gray and Hamann 2012),
subspecies approach (e.g. Pearman et al. 2010). Results of although they rely on the assumption that community
our approaches including intraspecific information concur composition and thereby ecological structure will endure.
closely with predictions of net change in the extent of Gray and Hamann (2011) found that if P. contorta popu-
ª2013TheAuthors.PublishedbyBlackwellPublishingLtd. 443
ImportanceofIntraspecificVariationinSDM B.Oneyetal.
120 (a) datasets that are more expansive and representative of the
species subject, and focus on the growth response rather
bitat (%) 100 tahnaanlomg otortaPl.itcyo.nItnortthae, tchaeseraonfgPeinoufstshyelvsepsterciise,sththeaEtuBreansiiatno
Ha Garzo´n et al. (2011) sampled represents a relatively small
e 80 portion of species range, and the performance of their
bl
uita conventional SDM may be a result of the small geo-
S 60 graphic range of the occurrence data; Thuiller et al.
e
g (2004) present results demonstrating that the smaller the
a
ent 40 range, the more likely the outlook produced by SDMs will
c MaxEnt 1 species
Per MaxEnt 3 subspecies be grim. Although tree mortality is directly related to spe-
Universal Transfer Function cies distributions, tree mortality is difficult to predict
20
(van Mantgem et al., 2009), as is also apparent in the
low-explained survival variance presented by Benito
120 (b)
Garzo´n et al. (2011), and influenced by a range of factors
%) inaddition to climate, that is, barkbeetles, over-browsing,
at ( 100 and other more proximal influences (Austin 2007). The
bit convergent results of the study by Benito Garzo´n et al.
a
H
e 80 (2011) and ours indicate the robustness of the underlying
bl phenomenon and highlight the utility of provenance trials
a
uit to investigate sub-species variation.
S 60
e RegardingtheuseofconventionalSDMs,themostcost-
g
nta and time-efficient method to incorporate intraspecific
e 40
c variability is to model sub-specific variants individually
er
P and combine the respective projections afterward, which
20 is presented first in Pearman et al. (2010) and again in
this study. This approach is of course less informative
1975 2020 2050 2080
than carrying out provenance trials, but serves the pur-
Climatic Period pose of creating more realistic species distributions than
conventional SDMs. Given the similarity of the methods
Figure4. PredictedsuitablehabitatareaofPinuscontortarelativeto
presented here that incorporate intraspecific variability
the reference climatic period 1950–2000 assuming (a) full dispersal
(subspecies model vs. UTF), it appears that “sub-clade”
and (b) no dispersal for each climatic period for all modeling
techniques. Modeling techniques which incorporate intraspecific (Pearman et al. 2010) or subspecies models are a more
variability – the UTF and subspecies model, predict more optimistic robust approach to using occurrence data. As Pearman
outcomesforPinuscontorta. et al. (2010) noted, the less prevalent conspecifics are
improperly represented in species (clade) models. Simi-
lationswereunabletomigrate,suitablehabitatwithinseed larly, although not quantitated, we find that the species
zones would decline by 14% (2020s), 22% (2050s) or 34% model underpredicts distributions of subspecies murraya-
(2080s), levels comparable to the “no dispersal” scenario na in the Sierra Nevada as well as subspecies contorta
of our intraspecific models (Fig. 4b). The species models along the coast, where they currently thrive, which is not
from McKenney et al. (2007) predict less habitat change the case in the subspecies model (Fig. 5). The niche of
than our species model, whereas Coops and Waring the most prevalent subspecies latifolia appears to be best
(2010) using a process-based model predict a similar represented in the species model. Pearman et al. (2010)
change in habitat change; nonetheless, the approaches uti- described prerequisites for dividing a species or taxon
lizing intraspecific information presented here predict into subcomponents, which summarize to ecogeographi-
much lesser climate change impacts, given both perfect cal distinction of within-taxon ecotypes, in other words,
and no migration. This further indicates that migration substantial niche differentiation which is spatially segre-
abilitiesshouldbeinvestigatedatsub-specieslevel. gated. Given the differences among the subspecies of
Benito Garzo´n et al. (2011) came to very similar con- P. contorta, the effects of climate may reasonably have a
clusions regarding the meaning of intraspecific informa- different effect on each subspecies and our results support
tion for Pinus sylvestris and Pinus pinaster in light of this notion (Fig. 5). Pearman et al. (2010) show that 7 of
climate change using provenance trial survival data. Our 10 species, which exhibit considerable intraspecific varia-
study complements and corroborates their results with tion, have greater projected range extents in future
444 ª2013TheAuthors.PublishedbyBlackwellPublishingLtd.
B.Oneyetal. ImportanceofIntraspecificVariationinSDM
120 (a) niches may also reflect evolutionary lag following the last
glaciation or recent climate change (O’Neill et al. 2008).
%) 100 Given the range of climates that P. contorta inhabits, it is
at ( very possible that the predicted suitable areas would host
bit P. contorta well.
a 80
H The work of Thomas et al. (2004) initiated a wide
e
bl discussion regarding the extinction risks associated with
a 60
uit global climate change. In reply, Harte et al. (2004) argued
S
e that the conclusions therein may be conservative in that
g 40
nta the averaging characteristic of SDM implies that one eco-
ce type is as capable as another ecotype, from a climatically
er 20 contorta
P different portion of the species’ range, of filling the
murrayana
latifolia modeledniche.ResultsfromtheIllingworthtrial(Rehfeldt
0
et al. 1999; Rehfeldt and Wykoff 2001; O’Neill et al. 2008)
120 (b)
and other large provenance trials (Reich and Oleksyn
2008)concur in that tree growth decreases with increasing
%) 100
bitat ( cHliomwaetvicerd,iwsthanencecboentswideeenrintgestthseiteUaTnFd(pi.oep.,uOlat’Nioenilloreitgianl..
a 80 2008) across all provenances, the contribution of intraspe-
H
e cific variation appears to buffer climate change impacts.
bl
uita 60 Thuilleret al.(2004)presentresultssuggestingthatspecies
S exhibiting greater niche variation have a greater chance of
e
g 40 copingwithclimatechange.Furthermore,resultsfromthis
a
nt study (as well as Pearman et al. 2010) indicate that con-
e
c
er 20 ventional species distribution modeling techniques are
P
unable to properly estimate the niche breadth of intraspe-
0 cific variants, especially if the intraspecific variants are
differently prevalent, and may therefore result in future
1975 2020 2050 2080 projections that overestimate suitability distribution losses
Climatic Period (e.g., Thomas et al. 2004; Thuiller et al. 2008; Urban et al.
2012). Considerable intraspecific variability is a common
Figure5. Relative predicted suitable habitat area of Pinus contorta
feature of plant species with broad distributions as indi-
subspecies assuming (a) full dispersal and (b) no dispersal for each
cated by local adaptations (e.g., Joshi et al. 2001; Hufford
climaticperiod(takenfromthesubspeciesmodel).
and Mazer 2003; Savolainen et al. 2007; Kreyling et al.
climates. Furthermore, their models, which incorporate 2012b). We find that the individual subspecies models
within-species variation, outperform their species models, show the same or better performance than the single
similar to our results (Fig. 2 and see Supporting Informa- species model (Fig. 4), in accordance to the results of
tion). Species with much intraspecific variability such as Kadmon et al. (2003). However, the performance of the
P. contorta exhibit niche diversity and breadth, which is combined subspecies model cannot easily be explained by
difficult to capture in a single model. In the case of the often better performance of models for species with
P. contorta, it appears that the niche breadth was underes- smallerniches,becausethecombinedmodeldoesnotnec-
timated by the species model. All models predict suitabil- essarily represent small niches, but rather a wider range of
ity in Eastern Canada and Southwest Alaska, where niches, and it is therefore relevant that the overall perfor-
P. contorta occurs rarely or is absent, suggesting that sub- mance of the combined model was not lower than that of
stantial niche overlap with other boreal tree species. As the single species model. Species with restricted distribu-
thoroughly discussed in the field of SDM (e.g., Guisan tions are on average more vulnerable to climate change
and Zimmermann 2000), only the realized niche is cap- (Thuiller et al. 2004; Schwartz et al. 2006), and species
tured by SDM approaches. Perhaps P. contorta could with littleintraspecific variability or niche breadth arealso
occur more frequently in many other parts of Canada mostlikelytobedisadvantaged.Whilethisstudyaddressed
and in southern Alaska, but natural history, processes intraspecific variation, all approaches assumed that the
such as dispersal and competition, and/or forestry man- derived relationships to the environment do not evolve, as
agement have shaped current observation patterns. theylikelywouldforrealpopulations(BenitoGarzo´net al.
Furthermore, difference in realized and fundamental 2011). However, our findings indicate that intraspecific
ª2013TheAuthors.PublishedbyBlackwellPublishingLtd. 445
ImportanceofIntraspecificVariationinSDM B.Oneyetal.
variation alone already allows for buffering against envi- Austin,M.2007.Speciesdistributionmodelsandecological
ronmental change. This is relevant as examples of niche theory:acriticalassessmentandsomepossiblenew
conservatism,thatis,thestaticnatureofaniche,occurfre- approaches.Ecol.Model.200:1–19.
quently in the fossil record (Davis and Shaw 2001) and in Beierkuhnlein,C.,D.Thiel,A.Jentsch,E.Willner,and
currentstudies(Wienset al.2009).Arguably,O’Neillet al. J.Kreyling.2011.EcotypesofEuropeangrassspeciesrespond
(2008) provide a compelling example of niche conserva- differentlytowarmingandextremedrought.J.Ecol.
tism in P. contorta: the best predictor of transplanted tree 99:703–713.
growth was the mean cold monthly temperature of its BenitoGarzo´n,M.,R.Alía,T.M.Robson,andM.A.Zavala.
provenanceclimate. 2011.Intra-specificvariabilityandplasticityinfluence
In conclusion, this study emphasizes the need to potentialtreespeciesdistributionsunderclimatechange.
investigate intraspecific variation, if applicable, when Glob.Ecol.Biogeogr.20:766–778.
considering assisted migration (Aitken et al. 2008; Best,A.S.,K.Johst,T.Münkemüller,andJ.M.J.Travis.
Kreyling et al. 2011). Whether using provenance trial data 2007.Whichspecieswillsuccessfullytrackclimatechange?
(e.g. O’Neill et al. 2008; Benito Garzo´n et al. 2011) or Theinfluenceofintraspecificcompetition anddensity
dependentdispersalonrangeshiftingdynamics.Oikos
acknowledging ecogeographically unique intraspecific
116:1531–1539.
variants with a conventional SDM technique (e.g., Pear-
Clark,J.S.,D.M.Bell,M.H.Hersh,andL.Nichols.2011.
man et al. 2010), intraspecific information is a necessary
Climatechangevulnerabilityofforestbiodiversity: climate
addition to sapient SDM analyses. Our results show that
andcompetition trackingofdemographicrates.Glob.
incorporation of intraspecific variation results in very dif-
ChangeBiol.17:1834–1849.
ferent, in this case much more positive, projections of
Coops,N.C.,andR.H.Waring.2010.Aprocess-based
future suitable habitat for a highly variable species.
approachtoestimatelodgepolepine(PinuscontortaDougl.)
Assessing the extent to which species are able to utilize
distributioninthePacificNorthwestunderclimatechange.
intraspecific variation presents a pertinent research chal-
ClimaticChange105:313–328.
lenge.
Davis,M.,andR.Shaw.2001.Rangeshiftsandadaptiveresponses
toQuaternaryclimatechange.Science292:673–679.
Acknowledgments Davis,M.,R.Shaw,andJ.Etterson.2005.Evolutionary
responsestochangingclimate.Ecology86:1704–1714.
The P. contorta distribution data in BC was provided by
Dobrowski,S.2011.Aclimaticbasisformicrorefugia:theinfluence
Nicholas Ukrainetz of the BC Forest Service. The authors ofterrainonclimate.Glob.ChangeBiol.17:1022–1035.
thank the US Forest Service Forest Inventory and Analysis
Dormann,C.F.2007.Promisingthefuture?Globalchange
Program for providing US P. contorta occurrence data.
projectionsofspeciesdistributions.BasicAppl.Ecol.
Numerous online resources were gratefully used to sup- 8:387–397.
plement occurrence data (See Supporting Information for Dormann,C.F.,O.Purschke,J.R.Márquez,S.Lautenbach,
thorough list). BR and JK acknowledge support by the andB.Schröder.2008.Componentsofuncertaintyin
“Bavarian Climate Programme 2020” within the joint
speciesdistributionanalysis:acasestudyofthegreatgrey
research center FORKAST. This publication was funded shrike.Ecology89:3371–3386.
by the German Research Foundation (DFG) and the Uni- Elith,J.,S.J.Phillips,T.Hastie,M.Dud(cid:1)ık,Y.E.Chee,andC.J.
versity of Bayreuth in the funding programme Open Yates.2011.AstatisticalexplanationofMaxEntforecologists.
Access Publishing. Divers.Distrib.6:1–15.
Fazekas,A.J.,andF.C.Yeh.2006.Postglacialcolonization
andpopulationgeneticrelationshipsinthePinuscontorta
Conflict of Interest
complex.Can.J.Bot.84:223–234.
None declared. Fielding,A.H.,andJ.F.Bell.1997.Areviewofmethodsfor
theassessmentofpredictionerrorsinconservation
presence/absencemodels.Environ.Conserv.24:38–49.
References
Fischer, D., S. M. Thomas, F. Niemitz, B. Reineking, and
Aitken,S.N.,S.Yeaman,J.A.Holliday,T.Wang,and C. Beierkuhnlein. 2011. Projection of climatic suitability
S.Curtis-McLane.2008.Adaptation,migrationor for Aedes albopictus Skuse (Culicidae) in Europe
extirpation:climatechangeoutcomesfortreepopulations. under climate change conditions. Global Planet. Change
Evol.Appl.1:95–111. 78:54–64.
Araju´o,M.B.,andM.Luoto.2007.Theimportanceofbiotic GRASSDevelopmentTeam.2012.GeographicResources
interactionsformodellingspeciesdistributionsunder AnalysisSupportSystem(GRASSGIS)Software.Open
climatechange.Glob.Ecol.Biogeogr.16:743–753. SourceGeospatialFoundation,USA.
446 ª2013TheAuthors.PublishedbyBlackwellPublishingLtd.