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

Learn More →

Reconstructing range dynamics and range fragmentation of European bison for the last 8000 years

Reconstructing range dynamics and range fragmentation of European bison for the last 8000 years Introduction Globally, biodiversity is in decline and conservationists face the challenge to identify pathways to preserving species and populations in an increasingly dynamic world ( Pressey , 2007 ; Brook , 2008 ; Franklin, 2010 ). The drivers of environmental change such as climate or land use change increasingly affect large regions at once. Designing resilient conservation strategies and deciding upon the most effective use of limited conservation funds therefore require moving beyond assessing individual sites towards range‐wide conservation planning at regional to continental scales ( Viña , 2010 ; Kuemmerle , 2011b ; Redford , 2011 ; Wikramanayake , 2011 ). Because many species often only persist in small, fragmented populations, this first and foremost necessitates to understand what constitutes species’ entire ranges and thus where species have occurred prior to large‐scale human disturbance ( Willis & Birks, 2006 ; Nogues‐Bravo, 2009 ). Moreover, assessing how species’ ranges have varied in the past can help conservation planners to better understand how species may respond to future environmental change ( Wiens , 2009 ; Graham , 2010 ). Climate and land use change are both rapidly accelerating, yet neither their relative importance nor their synergistic effects on biodiversity are well understood ( Brook , 2008 ). The increasing availability of climate, vegetation and land use reconstructions ( Kaplan , 2009 ; Schurgers , 2009 ; Klein Goldewijk , 2010a ) provides new opportunities to study past range dynamics and to disentangle the combined effects of climate and land use change on species’ distributions. Furthermore, assessing how expanding human populations and land use has diminished and fragmented species’ ranges can help conservationists understand the sensitivity of species to human disturbance, to identify critical tipping points in habitat fragmentation and to locate refugia ( Nogues‐Bravo, 2009 ; Graham , 2010 ; Swift & Hannon, 2010 ), all of which is important for designing resilient conservation strategies. Species distribution models (SDM) have become increasingly popular for mapping current or potential future ranges of species. SDM relate occurrence locations to a suite of environmental variables ( Elith & Leathwick, 2009 ) and thus assess species’ realized niches (i.e. the environmental conditions used by a species, Soberon & Nakamura, 2009 ). This can be problematic if a given species occupies only a portion of its potential range, for example because of past habitat loss, overhunting or dispersal barriers ( Nogues‐Bravo, 2009 ). Parameterizing SDM based on historic occurrence data prior to widespread human disturbance results in more robust assessments of species’ distributions in such cases ( Nogues‐Bravo, 2009 ). However, although such hindcasting approaches are powerful tools to better understand species’ response to environmental change and thus ultimately to improve conservation strategies ( Nogues‐Bravo, 2009 ; Willis , 2009 ; Graham , 2010 ), SDM hindcasting applications are still scarce. Large carnivores and herbivores are archetypical examples of species whose current distributions are strongly determined by human disturbance ( Mladenoff , 1999 ; Kuemmerle , 2010 ). These species require large tracts of intact habitat, often conflict with land use, and are prone to poaching ( Breitenmoser, 1998 ; Gordon & Loison, 2009 ). As a result, large carnivores and herbivores are among the most challenging species to preserve in human‐dominated landscapes, and their populations have declined world‐wide. Safeguarding existing populations of large herbivores and carnivores and restoring their crucial ecological roles are thus important and well‐accepted conservation goals that may also benefit many other species ( Carroll , 2001 ; Vera , 2006 ; Gordon & Loison, 2009 ). Europe’s temperate zone has a long history of human settlement and land use ( Bramanti , 2009 ; Kaplan , 2009 ), and this has resulted in the extirpation of large carnivores and herbivores throughout the majority of their former ranges ( Breitenmoser, 1998 ; Vera , 2006 ). The European bison (or wisent, Bison bonasus L.) is Europe’s largest land mammal and the last surviving large grazer, as well as a prime example of a species that today only persists in a small, fragmented population. Bison disappeared from the wild in the early 20th century, but a systematic breeding and reintroduction programme prevented their extinction. Today, wild European bison occur in about 35 small, isolated herds in Central and Eastern Europe ( Pucek , 2004 ; Krasinska & Krasinski, 2007 ). As genetic diversity of the European bison population is low, the species continues to be at risk ( Tokarska , 2011 ), and the survival of wild bison will depend on establishing large metapopulations ( Perzanowski , 2004 ; Kuemmerle , 2011a ). However, identifying appropriate sites to establish large bison metapopulations is challenging because of two reasons. First, the habitat requirements of European bison are not fully understood. European bison are traditionally thought of as a species associated with closed temperate forests ( Heptner , 1961 ), but such ecosystems may simply have been their last refuges. Likewise, the species is generally considered a grazer, although browsing could contribute substantially to their diet ( Kowalczyk , 2011 ). Second, substantial uncertainty exists regarding both the species’ historic distribution ( Heptner , 1961 ; Pucek , 2004 ) and the effects of past environmental change and expanding human populations on the bison’s range. Recently, human population density and land use intensity have declined in many regions across Central and Eastern Europe in the wake of the collapse of socialism ( Ioffe & Nefedova, 2004 ; Baumann , 2011 ). This may offer a window of opportunity to implement a broad‐scale conservation strategy for safeguarding European bison and other large carnivores and herbivores ( Kuemmerle , 2011b ). A better understanding of habitat characteristics and historic distributions of European bison is thus urgently needed to identify priority sites for European bison conservation. Here, our goal was to reconstruct European bison range dynamics for the last 8000 years. To parameterize our SDM, we used a comprehensive set of Holocene European bison occurrences from before the species went extinct in the wild together with bioclimatic variables and vegetation reconstructions from a dynamic vegetation model as predictors. Including vegetation reconstructions addresses two limitations of SDM based solely on bioclimatic variables: First, bioclimatic variables tend to result in simplistic niche characterizations and may overpredict geographic distributions substantially ( Pearson & Dawson, 2003 ; Wiens , 2009 ). Second, climate‐based assessments of range shifts for species at higher trophic levels disregards that changes in vegetation communities, on which these species may depend, are slower than changes in climate. Dynamic vegetation models can capture such vegetation transitions ( Schurgers , 2009 ; Hickler , 2011 ). Using simulated vegetation in tandem with niche models can therefore improve range reconstructions for species at higher trophic levels, but to our knowledge, no prior study has done this. Specifically, we asked the following research questions: 1 What characterized the realized niche of European bison during the Holocene? 2 What was the geographic distribution of European bison during the last 8000 years and how did this distribution vary? 3 How have human population growth and farmland expansion affected the habitat of European bison? Methods Holocene records of European bison occurrence We used two types of historic European bison occurrences: (1) bone findings from archaeological sites and natural deposits, and (2) written records of bison occurrences. Locations of bison bones stemmed from a comprehensive archaeozoological database containing more than 7000 faunal assemblages from 31 countries across Europe since the last glaciation ( Benecke, 1999 ). This database contains 169 geolocated and dated findings of European bison remains ( Benecke, 2005 ), mainly from the mid‐ and late Holocene ( Fig. 1 ). Because the Benecke (1999) database is particularly rich in faunal assemblages from Southern, Western, Northern and Central Europe, but relatively sparse in Eastern Europe, we also included 45 occurrence locations from Heptner (1961) covering the European region of the former Soviet Union ( Fig. 1 ). These occurrence locations were derived from written records from various historical sources from the Middle Ages until the 20th century ( Heptner , 1961 ). Following Benecke (2005) , we assigned all occurrence locations to one of three time periods: (1) mid‐Holocene (6000 bc –1000 bc ), (2) Iron Age (1000 bc –600 ad ) and (3) Middle Ages and Modern Period (600 ad –1900 ad ). 1 European bison occurrences during the mid‐ and late Holocene based on bone findings ( Benecke, 2005 ) and written records ( Heptner , 1961 ). Only bison occurrences before the extirpation of the species in the wild in the early 20th century were used. The archaeozoological database ( Benecke, 1999 ) contains numerous faunal assemblages from Western and Southern Europe, but bison remains were absent from these regions during the mid‐ and late Holocene. (Projection: Albers equal‐area conic). Predictor variables Our study region covered all of Europe (15° W – 65° E and 75° N – 30° N). Climate variables for the dynamic vegetation model and the range reconstructions were derived from a transient simulation from 7000 bc to 2000 ad ( Schurgers , 2006 ) with a complex earth system model that includes atmosphere and ocean dynamics, as well as terrestrial and marine carbon cycling ( Mikolajewicz , 2007 ). The anomalies (averaged per 100‐year period) from this simulation compared to a 10,000‐year control simulation were superimposed on a detrended climate dataset with monthly average temperature, precipitation and cloud cover for the 20th century ( Mitchell & Jones, 2005 ). We then derived annual precipitation, potential evapotranspiration, growing degree days (temperature sum above a base temperature of 5°), as well as the minimum, mean and maximum temperature of the last year of each 100‐year period. Winter severity is a critical factor affecting ungulate survival, including that of European bison ( Krasinska & Krasinski, 2007 ). As a proxy variable for winter severity, we calculated the mean temperature of the coldest quarter (December–February). Vegetation reconstructions were simulated using an updated version of the dynamic vegetation model LPJ‐GUESS (v2.1, Smith , 2001 ). Vegetation was represented by 11 global plant functional types (PFTs) with different bioclimatic niches, growth forms, leaf phenology type (i.e. evergreen, summergreen or raingreen) and photosynthetic pathway (C3 or C4). The LPJ‐GUESS model successfully reproduces the main features of the global distribution of vegetation types, in particular in the study area. In this study, the model was run continuously for 10,000 years, using the first 1000 years to spin up the vegetation from bare ground. Model inputs consisted of the climate data described above, present‐day soil texture ( Sitch , 2003 ), and reconstructed and observed CO 2 concentrations ( Indermühle , 1999 ). As input variables for the SDM, we extracted the leaf area index (LAI; averaged per 100‐year period) for the PFTs that occurred in our study region. This yielded seven vegetation variables: LAI of boreal coniferous trees, LAI of temperate, broadleaved deciduous trees, LAI of temperate, shade‐intolerant deciduous trees, LAI of temperate evergreen trees, LAI of C3 grasses, total LAI and total woody LAI (see Fig. S1 in Supporting Information). All climate and vegetation variables were projected to the Albers equal‐area conic coordinate system. Reconstructing European bison range dynamics To characterize the realized niche of European bison during the Holocene, we used multitemporal niche calibration ( Nogues‐Bravo, 2009 ). This approach accounts for the fact that the environmental conditions a species utilizes may vary through time by matching each occurrence point with the environmental data from the time when a species occurred at a location. Using the resulting multitemporal dataset, we parameterized a single niche model that was then projected onto the environmental layers of each time period. Thus, all information about environmental conditions a species utilized throughout the full time span considered (6000 bc –2000 ad ) was used to describe its geographic distribution at a particular point in time ( Nogues‐Bravo, 2009 ). Calibrating niche models in this way results in more accurate and robust range reconstructions than reconstructions based on either past or present environmental predictors alone ( Broennimann & Guisan, 2008 ; Nogues‐Bravo, 2009 ). We extracted environmental conditions for the bison occurrences using the climate and vegetation layers closest to the mid‐point of our time periods (3000 bc , 0 ad , 1000 ad for the mid‐Holocene, Iron Age, and Middle Ages, respectively). We used maximum entropy modelling ( Phillips , 2006 ), to reconstruct the Holocene range of European bison. The maximum entropy approach assumes that the potential geographic distribution of a species is a probability distribution π over all locations (i.e. cells). To approximate π, a probability distribution is derived, considering constraints imposed by the environmental predictors at the occurrence locations ( Elith , 2011 ). The distribution with maximum entropy approximates π best because it is least constrained. Maximum entropy modelling requires only presence data, works well with small sample sizes and is among the highest performing SDM approaches ( Elith , 2006 ). A detailed mathematical description of maximum entropy modelling is provided in Phillips (2006) and Elith (2011) . To fit maximum entropy models, we used Maxent (v3.3.3). All model runs used default regularization, a maximum of 2500 iterations and 10,000 random background points ( Phillips & Dudik, 2008 ). Background points were randomly sampled from the environmental layers from all three time periods of our occurrence data. To prevent overfitting, we used only quadratic and hinge features. We built two groups of Maxent models: (1) models based on climate variables only and (2) models based on climate and vegetation variables. Furthermore, a few of our predictor variables were collinear ( r > 0.7) and we fitted alternative models in such cases, retaining the variable yielding higher goodness‐of‐fit. We assessed the goodness‐of‐fit of our models based on the receiver operating characteristics (ROC), which compares the true positive rate (i.e. sensitivity) versus the false positive rates (i.e. 1−specificity) across all possible thresholds, and by calculating the area under the curve (AUC), which provides a threshold‐independent measure of model performance ( Phillips , 2006 ). We used a fivefold cross‐validation strategy (i.e. five model runs based on 80% of the occurrence data, while retaining 20% for validation) and calculated mean AUC values and AUC standard errors for each model run. Final range maps were calculated as the average of the five replicate runs per millennia, and we used a logistic link function to yield a relative environmental suitability index (ESI) between zero and one ( Phillips & Dudik, 2008 ). We also used a jackknife procedure to assess variable importance, where we fitted models excluding a particular variable as well as single‐variable models and compared changes in AUC compared to the full model. Once Maxent models for both model groups were parameterized, we projected bison ranges for each millennia from 6000 bc to 2000 ad . Because European bison are a non‐migratory species ( Krasinska & Krasinski, 2007 ), our range maps depict the potential year‐round occurrence of European bison. To summarize the range maps, we extracted the area covered by ESI values larger than the minimum and the 10‐percent presence ESI values at our occurrence (i.e. training) locations for each millennium, and we calculated the median‐filtered average ESI map across all millennia. Reconstructing habitat fragmentation We used the HYDE 3.1 database that provides human population and land use grids for the last 12,000 years at a spatial resolution of 5 arc‐min ( Klein Goldewijk , 2010a,b ). We calculated the relative share of cropland and pasture for each grid cell and used human population density for all years for which ESI maps had been calculated (every millennium since 6000 bc ) plus the years 1800 ad and 1900 ad . All maps were transformed to the Albers equal‐area conic projection and resampled to 10‐km resolution. To assess how population and land use expansion affected the bison’s range, we crosstabulated the cropland, pasture and population density maps with the ESI maps. We used 5%‐wide bins for the cropland and grassland maps, logarithmic bins for the population density map (using the breaks 1, 10, 100, 1000, 10,000 and 100,000), and the minimum and 10% presence thresholds for the ESI maps. For 1800 ad and 1900 ad , we used the 2000 ad ESI map. To assess how human population and land use expansion fragmented the bison’s range, we masked all areas above the cropland, pasture and human population thresholds that were assumed to exclude bison populations. We compared eight scenarios: two levels of ESI thresholds (minimum and 10% presence) and four levels of land use and population thresholds: (1) cropland or pasture density > 5%/population density > 5 persons km −2 , (2) > 10%/> 10 persons km −2 , (3) > 15%/> 20 persons km −2 , and (4) > 20%/> 40 persons km −2 . Next, we derived the remaining habitat area, the number of habitat patches, mean patch size, the largest patch index (i.e. percentage of the landscape comprised by the largest habitat patch), and the mean proximity index, which measures isolation of patches ( McGarigal , 2002 ). As a neighbourhood distance for the proximity index, we used 300 km (i.e. roughly the maximum recorded dispersal distance of European bison, Krasinska & Krasinski, 2007 ; Kuemmerle , 2011a ). Results European bison inhabited a broad range of environmental conditions during the last 8000 years ( Fig. 2 ). Areas inhabited by bison had a range of average annual temperature between 0 and 13.2 °C, and bison tolerated average winter (i.e. December–February) temperatures as low as −14.9 °C. Bison also tolerated a wide range of annual precipitation, but avoided dryer regions (e.g. steppe regions). Bison occurrences were characterized by high vegetation cover and a mixture of woody vegetation and grasses. We did not find a clear preference for broadleaved versus coniferous forests, but bison occurred mainly in forests consisting of shade‐tolerant species ( Fig. 2 ). 2 Range of climate and vegetation conditions covered by bison occurrence locations during the mid‐ and late Holocene. Climate and vegetation information for the individual bison occurrence points (dated bone findings or written records) were extracted from the time periods when bison occurred at these locations and pooled to allow for a comprehensive description of the bison’s realized niche for the time period 6000 bp –2000 ad . Our final Maxent model induced three climate variables (average annual temperature, average winter temperature, and annual precipitation) and six vegetation variables (LAIs of boreal coniferous trees, temperate broadleaved trees, temperate shade‐intolerant broadleaved trees, temperate evergreen trees, C3 grasses and total LAI). The goodness‐of‐fit of this model was high, with an average cross‐validated AUC value of 0.920 and a standard error of 0.012. Based on our jackknife analysis, LAI of temperate broadleaved and LAI of shade‐intolerant broadleaved trees were the most important variables in our model, followed by mean winter temperature (which also contained the most information not present in any other variable). Variable response curves confirmed that environmental suitability for European bison was high in areas of high total vegetation cover that were dominated by temperate broadleaved forests but also contained grasses. Our models improved when including vegetation predictors compared to using climate predictors only. Our best model based on climate predictors alone included average annual temperature, average winter temperature and annual precipitation (AUC = 0.901). Distributional patterns differed substantially between the two model types, generally resulting in fewer and more clustered areas of high ESI values when including also vegetation predictors. For example, extensive areas in southern Scandinavia and Turkey appeared climatically suitable for European bison, but had low ESI values because of their vegetation composition (see Fig. S2). Projecting our calibrated niche model to all nine millennia and calculating average ESI values revealed that the heartland of European bison during the mid‐ and late Holocene was in Central and Eastern Europe with highest suitability values in contemporary Ukraine, Poland, Belarus, the Baltics, Romania and European Russia. ESI values were also moderately high in southern Scandinavia, the Balkans, parts of Turkey, and the Caucasus. In contrast, Western Europe contained few suitable areas ( Fig. 3 ). 3 Distribution of European bison during the mid‐ and late Holocene [average environmental suitability index (ESI) from 6000 bc –2000 ad ]. ESI values were derived by projecting a time‐calibrated niche model based on Holocene bison occurrences and a maximum entropy model to the environmental layers (climate and vegetation reconstructions) for each millennium. (Projection: Albers equal‐area conic). The geographic distribution of suitable areas was relatively stable during the last 8000 years ( Fig. 4 ). Variability in time occurred mainly in European Russia and southern Scandinavia, with a more northern distribution before 1000 bc . Today’s distributional patters are most similar to those from the mid‐Holocene (6000 bc and 5000 bc ). The distribution of high ESI values was not contiguous, but limited to isolated patches occurring, for example, in southern France and the Caucasus ( Fig. 4 ). Some of these patches, for instance the Caucasus, were not always isolated though, and the Caucasus did harbour a subspecies of European bison until the early 20th century ( Krasinska & Krasinski, 2007 ). The area covered by ESI values above the minimum and 10% presence thresholds did also not vary substantially (2.63 and 5.85 million km 2 on average, respectively, Fig. 5 ). 4 Dynamics in the distribution of environmentally areas suitable for European bison during the Holocene based on projecting a time‐calibrated niche model to the environmental layers of each time period. (Projection: Albers equal‐area conic). 5 Changes in the area characterized by environmental suitability index (ESI) values larger than the minimum value and the 10%‐percentile value observed at any of the locations where European bison were present during the mid‐ and late Holocene. Comparing our distributional maps with maps of cropland, pasture, and human population density suggests that the area suitable and available for European bison diminished steadily during the last 8000 years. For example, of the entire area with ESI values above the minimum presence threshold, almost 60% contained cropland ( Fig. 6a ) and more than 70% had a population density exceeding 10 persons km −2 ( Fig. 6e ) by 0 ad . Cropland expansion affected the bison’s range more than pasture expansion, although most areas with high ESI values contained pasture land by 1000 ad ( Fig. 6c,d ). After 1800 ad , areas with high ESI and low cropland, pasture and human population density decreased dramatically. For instance, of all the areas with ESI values above the 10% presence threshold, only 1.5% were without farmland in 2000 ad ( Fig. 6b ) and only 0.22% had human population densities < 10 persons km −2 ( Fig. 6f ). 6 Distribution of environmentally suitable areas for European bison among areas with varying cropland (a and b), pasture (c and d) and population densities (e and f). Left column: 100% = all areas with environmental suitability index (ESI) values larger than the minimum presence threshold; right column: 100% = all areas with ESI larger than the 10%‐presence threshold. Land use and human population expansion gradually fragmented the area suitable for European bison until about 1000 ad , after which fragmentation increased sharply ( Fig. 7 ). For example, the number of habitat patches increased steadily between 6000 bc and 1000 ad ( Fig. 7a ), while mean patch area ( Fig. 7b ) and largest patch index ( Fig. 7c ) decreased, indicating fragmentation. Range fragmentation reached a threshold in 1000 ad , whereafter all fragmentation parameters changed drastically. Our sensitivity analyses revealed that this general pattern did not depend on a particular farmland or human population density threshold (see Fig. S3). Bison habitat was more fragmented and less widespread when applying the 10% presence compared to the minimum presence ESI threshold ( Fig. 7 ). 7 Number of habitat patches, mean patch area, largest patch index and mean proximity index of the remaining European bison habitat after masking areas with high cropland, pasture or human population density. We masked areas with cropland or pasture density > 10% or human population density > 10 persons km −2 and used two ESI thresholds to derive habitat for each time layer: minimum presence (dashed line) and 10% presence (straight line) ESI values. Discussion Using a time‐calibrated niche model in tandem with a dynamic vegetation model, we found that European bison habitat preferences may be broader than previously assumed. During the mid‐ and late Holocene, the heartland of European bison was in Central and Eastern Europe, and our models suggested that bison had a more eastern and northern distribution than assumed by expert‐based assessments. Whereas climate and vegetation transformations did not affect the distribution of European bison substantially during the last 8000 years, human population growth and the expansion of farming drastically diminished and fragmented European bison habitat, likely reaching a tipping point during the last 1000 years. Our findings have important conservation implications, suggesting that efforts to preserve bison should move beyond the paradigm that bison are a species of broadleaved forests only and focus on the northern and eastern edge of its current distribution to establish large bison metapopulations. More generally, our study showed that hindcasting approaches can help to better characterize species niches and to reconstruct range dynamics, which is crucial to design effective conservation strategies that cover the entire ranges of endangered species. Traditionally, bison were considered to thrive only in interior temperate broadleaved forests, but this may simply be because such habitats were the species’ last refuges ( Pucek , 2004 ; Krasinska & Krasinski, 2007 ; Kerley , 2011 ). Our analyses showed that bison may utilize a wide range of forest types and occupy temperate broadleaved and southern boreal mixed forests alike. These findings are in line with both field studies ( Perzanowski , 2008 ), fine‐scale habitat assessments ( Kuemmerle , 2010, 2011b ) and expert‐based assessments ( Benecke, 2005 ). Likewise, the diet preference of European bison is not fully understood. Previous studies found browsing to play a minor role ( Gebczynska , 1991 ; Krasinska & Krasinski, 2007 ), and physiological studies suggest that bison are mainly grazers ( Mendoza & Palmqvist, 2008 ). Yet, our models showed that bison selected for areas containing both forest and grass cover, but were most widespread during the mid‐ and late Holocene in closed forests (as opposed to more open forests, Fig. 2 ). This suggests that browsing may have accounted for a substantial share of bison’s diet, confirming a recent field study finding a browsing share of up to 65% if supplementary feeding is excluded ( Kowalczyk , 2011 ). While it was interesting to see that we did not find support for bison inhabiting grass‐dominated regions (e.g. steppes in contemporary southern Russia and parts of Ukraine) during the mid‐ and late Holocene, we would like to point out that we assessed only the bison’s realized niche. The species’ fundamental niche could be wider, possibly including grassland conditions, especially when considering that European bison evolved from the steppe bison ( Bison priscus ) ( Pucek , 2004 ; Benecke, 2005 ). Likewise, although vegetation dynamics during the mid‐ and late Holocene did not affect bison distribution substantially, the transformation from grasslands to forest‐dominated ecosystems during the early Holocene likely affected bison distribution ( Pucek , 2004 ; Benecke, 2005 ). The Holocene distribution of European bison that our model predicted differed substantially from previous, expert‐based assessments of the bison’s range ( Heptner , 1961 ; Pucek , 2004 ; Benecke, 2005 ; Sipko, 2009 ; Tokarska , 2011 ). We found the heartland of European bison to be in Central and Eastern Europe (i.e. consistent with prior assessments), but the range of European bison varied strongly on its eastern and northern edge during the last 8000 years. These are also the regions where expert‐based range assessments disagree strongest. Overall, our analyses suggested that suitable regions for European bison stretched further north and east than previously thought, at least during parts of the Holocene. Bison’s range in these regions was limited by winter severity, the main factor determining ungulates survival in general and the most important climate variable in our models, which supports the ecological realism of our model. We also note that there is evidence that bison may have indeed inhabited these more northern regions ( Sipko, 2009 ). European bison distribution during the Holocene did not extend substantially into Western Europe (e.g. France or northern Spain). While this contrasts with most prior range assessments (but see Tokarska , 2011 ), and although the species was present there during the Late Stone Age (e.g. as documented in the cave paintings of Altamira and Lascaux) and possibly also the early Holocene, no bison remains from the mid‐ and late Holocene have been found there despite a large number of archaeozoological assemblages ( Vigne , 2003 ; Benecke, 2005 ). Some of these regions were grass‐dominated (e.g. Spain), had a high share of evergreen broadleaved trees (e.g. southern France), or forests with very high LAIs (and low grass cover), all of which differed substantially from the environmental conditions at our occurrence locations ( Fig. 2 ). Overall, this suggests that European bison were rare in Western and South‐western Europe already by the mid‐Holocene. We caution, however, that the reasons for this remain somewhat unclear, as land use and human disturbance may also have pushed bison out of otherwise suitable areas on the south‐western and western edge of their range (see below). Potential causes of the extirpation of European bison from large parts of their range are changing environmental conditions, habitat fragmentation and overexploitation. Environmental change has triggered range contractions and population collapses of large mammals in Northern Eurasia during the Holocene ( Schmölcke & Zachos, 2005 ; Nogues‐Bravo , 2008 ), but our results confirm that human factors were the main cause of the decline of European bison. Areas suitable for European bison were relatively stable both in geographic extent and location. On the other hand, the expansion of settlements and farmland during the Holocene severely reduced and fragmented European bison habitat, particularly in Southern and Central Europe. The expansion of human populations and farming may also explain why some suitable regions within the bison’s range were apparently not occupied during the mid‐ and late Holocene, such as the Balkans, Bulgaria, southern Romania or southern France. Settlements and farming were already widespread in these regions during the mid‐Holocene ( Bramanti , 2009 ), and forest cover was likely drastically reduced by 1000 bc ( Kaplan , 2009 ). On the other hand, later farmland expansion and lower human density may be the reason why bison could hold out longest in the north of Central and Eastern Europe ( Heptner , 1961 ), a region that maintained high forest cover until the 19th century ( Kaplan , 2009 ). Habitat loss and fragmentation did not occur uniformly in time, though. Both the extent and the connectivity of habitat plummeted during the last 1000 years, particularly since 1800 ad ( Fig. 7 ), when habitat availability was well below thresholds that may result in extirpation ( Fahrig, 2002 ). This suggests that the European bison population as a whole may have crossed a tipping point towards population collapse. Habitat fragmentation was also exacerbated by increased hunting, and overexploitation, which undoubtedly contributed to the decline and eventual extirpation of European bison. Our models yielded high goodness‐of‐fit and plausible distributional patterns that are consistent with independent assessments of the range contraction and population decline of European bison during the Holocene. A few sources of uncertainty remain. First, while we used a large set of occurrence data from both bone findings and written records, we cannot fully rule out bias in our occurrence data (e.g. because bone preservation may differ among soil types, sampling effort may be unequally distributed, B. bonasus remains cannot be clearly identified or bison did not occupy all environmentally suitable regions, Schmölcke & Zachos, 2005 ). While faunal assemblages were also widespread in regions where no bison records were found (e.g. France, Spain), bone findings were often associated with human settlements, adding uncertainty in areas that were settled later ( Bramanti , 2009 ; Kaplan , 2009 ). However, sampling bias does not affect our range reconstructions as long as environmental analogues (in either time or space) to undersampled regions exist in our occurrence dataset, representing a major advantage of our time‐calibrated niche model over mono‐temporal SDM or traditional range reconstructions based on visual interpretation of historic occurrence locations. Second, we used only a single climate and vegetation reconstruction. However, climate variations over this period are primarily driven by well understood orbital changes, and the main determinants of vegetation distribution and structure in Europe are well known and covered by LPJ‐GUESS ( Hickler , 2011 ). We also note that our vegetation constructions compare favourably to pollen reconstructions ( Prentice , 1998 ; Tarasov , 1998 ) and that we are not aware of any alternative climate reconstructions from Global Climate Models or transient vegetation reconstructions for the time period we assessed. Third, while our cross‐validation strategy allows for robust comparisons among alternative models, it may be optimistic in assessing overall model fit ( Araújo , 2005 ). However, cross‐validation remains the only feasible approach to assess the goodness‐of‐fit of our models as an independent set of presence/absence data cannot be gathered retrospectively ( Araújo , 2005 ). Last, while our sensitivity analyses showed that our main conclusions remain unaffected by the particular choice of land use and human population thresholds, the HYDE database has some uncertainties, especially regarding pastures and for older time layers ( Klein Goldewijk , 2010a,b ). Comparisons with other land use reconstructions suggest relatively consistent farmland expansion and deforestation patterns ( Olofsson & Hickler, 2008 ; Kaplan , 2009 ; Gaillard , 2011 ). Yet, the HYDE scenario likely represents a conservative estimate of the magnitude of land use change ( Kaplan , 2009 ), meaning that range fragmentation could have occurred even earlier and more rapidly than indicated by our results. Our range results have a number of important implications for the conservation of European bison and other large herbivores and carnivores. First, our results showed that the habitat preferences of European bison during the Holocene were broader than previously thought, with bison thriving in semi‐open areas and in broadleaved, mixed and coniferous forests alike ( Fig. 2 ). This provides further support for views that bison conservation should also focus on areas outside deciduous broadleaved forests and that the area of potential suitable European bison habitat may be large ( Pucek , 2004 ; Kerley , 2011 ; Kuemmerle , 2011b ). Second, our hindcasting approach suggests that European bison had a larger range in Europe’s north and east than previously assumed, yet that bison were absent or only present in very low numbers in Western Europe already during the mid‐ and late Holocene. As preserving European bison in the wild will depend on establishing large metapopulation ( Pucek , 2004 ), efforts to establish such metapopulations should focus especially on central European Russia and the Baltic states. Rural populations and land use there have declined drastically there since 1991 ( Ioffe , 2004 ; Baumann , 2011 ), which could benefit the establishment of a large bison population. Moreover, whereas unmanaged grasslands are sparse in Western Europe, abandoned and fallow farmland is now widespread in Europe’s East, allowing bison to utilize the broad range of habitat types they have used in the past. Third, future reintroductions should also focus on the northern edge of the current European bison distribution. Our range reconstructions showed that bison extended substantially more northwards in the past, especially during warmer periods ( Fig. 4 ). Considering that climate change will shift northern European regions to more temperate conditions during the 21st century ( Hickler , 2011 ), northern European Russia may offer suitable candidate sites for establishing a large bison metapopulation ( Sipko, 2009 ; Kuemmerle , 2011b ). Fourth, our study showed that drastic habitat fragmentation likely contributed substantially to the decline of the European bison. As linking isolated bison populations is a conservation priority identified in the European bison species action plan ( Pucek , 2004 ), this provides another argument for focusing European bison conservation efforts on Belarus, European Russia and the Baltics, where habitat fragmentation is still much lower than in Western and Central Europe ( Angelstam , 1997 ). Fifth, conservation managers should interpret low ESI values in parts of Western Europe carefully, as there is a possibility that bison were pushed out of suitable habitat there prior to our assessment period. However, even if these regions contained suitable habitat in the past, we believe that the prospects for the large bison metapopulations needed to ensure the persistence of the species in the wild are low because of conflicts with land use ( Kuemmerle ., 2011b ). Finally, our assessments confirmed that during the last 8000 years, bison in the Caucasus were separated several times from the rest of the European bison range, suggesting that the two bison lines may represent different subspecies and that conservation managers should continue to keep these lines separated. Reconstructing range dynamics of European bison suggested that land use change was the major force contributing to the catastrophic decline of European bison, whereas environmental change did not affect the distribution of bison substantially during the last 8000 years. Most future range assessments currently focus exclusively on possible future climate effects. Our study thus emphasizes the need to consider land use change when assessing species’ ranges and, more generally, the future of biodiversity ( Brook , 2008 ; Leadley , 2010 ). While conservation should not be about mimicking the past, hindcasting can improve our understanding of species’ niches and historic distributions. Such baselines are urgently needed for range‐wide conservation planning and to help conservation managers to design resilient strategies to preserve species in an increasingly human‐dominated world. Acknowledgements We thank L. Baskin, W. Cramer, C. Fløjgaard, P. Daskiewicz, K. Perzanowski, T. Samojlik, T. Sipko, J.‐D. Vigne and D.M. Waller for valuable discussions, and N. Benecke and V.G. Heptner for compiling the comprehensive datasets of historic European bison occurrences. Associate editor N. Roura‐Pasqual, J.P.G.M. Cromsigt and an anonymous reviewer are thanked for constructive comments on prior manuscript versions. We gratefully acknowledge support by the Alexander von Humboldt Foundation, the EU project ECOCHANGE (FP6‐036866), the NASA Biodiversity and NASA Land Cover and Land Use Change programmes, and the LOEWE initiative for scientific and economic excellence of the German federal state of Hesse. This study is a contribution to the Swedish Strategic Research area Biodiversity and Ecosystem services in a Changing Climate (BECC). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Diversity and Distributions Wiley

Reconstructing range dynamics and range fragmentation of European bison for the last 8000 years

Loading next page...
 
/lp/wiley/reconstructing-range-dynamics-and-range-fragmentation-of-european-tXdt99ZEUQ

References (74)

Publisher
Wiley
Copyright
© 2011 Blackwell Publishing Ltd
ISSN
1366-9516
eISSN
1472-4642
DOI
10.1111/j.1472-4642.2011.00849.x
Publisher site
See Article on Publisher Site

Abstract

Introduction Globally, biodiversity is in decline and conservationists face the challenge to identify pathways to preserving species and populations in an increasingly dynamic world ( Pressey , 2007 ; Brook , 2008 ; Franklin, 2010 ). The drivers of environmental change such as climate or land use change increasingly affect large regions at once. Designing resilient conservation strategies and deciding upon the most effective use of limited conservation funds therefore require moving beyond assessing individual sites towards range‐wide conservation planning at regional to continental scales ( Viña , 2010 ; Kuemmerle , 2011b ; Redford , 2011 ; Wikramanayake , 2011 ). Because many species often only persist in small, fragmented populations, this first and foremost necessitates to understand what constitutes species’ entire ranges and thus where species have occurred prior to large‐scale human disturbance ( Willis & Birks, 2006 ; Nogues‐Bravo, 2009 ). Moreover, assessing how species’ ranges have varied in the past can help conservation planners to better understand how species may respond to future environmental change ( Wiens , 2009 ; Graham , 2010 ). Climate and land use change are both rapidly accelerating, yet neither their relative importance nor their synergistic effects on biodiversity are well understood ( Brook , 2008 ). The increasing availability of climate, vegetation and land use reconstructions ( Kaplan , 2009 ; Schurgers , 2009 ; Klein Goldewijk , 2010a ) provides new opportunities to study past range dynamics and to disentangle the combined effects of climate and land use change on species’ distributions. Furthermore, assessing how expanding human populations and land use has diminished and fragmented species’ ranges can help conservationists understand the sensitivity of species to human disturbance, to identify critical tipping points in habitat fragmentation and to locate refugia ( Nogues‐Bravo, 2009 ; Graham , 2010 ; Swift & Hannon, 2010 ), all of which is important for designing resilient conservation strategies. Species distribution models (SDM) have become increasingly popular for mapping current or potential future ranges of species. SDM relate occurrence locations to a suite of environmental variables ( Elith & Leathwick, 2009 ) and thus assess species’ realized niches (i.e. the environmental conditions used by a species, Soberon & Nakamura, 2009 ). This can be problematic if a given species occupies only a portion of its potential range, for example because of past habitat loss, overhunting or dispersal barriers ( Nogues‐Bravo, 2009 ). Parameterizing SDM based on historic occurrence data prior to widespread human disturbance results in more robust assessments of species’ distributions in such cases ( Nogues‐Bravo, 2009 ). However, although such hindcasting approaches are powerful tools to better understand species’ response to environmental change and thus ultimately to improve conservation strategies ( Nogues‐Bravo, 2009 ; Willis , 2009 ; Graham , 2010 ), SDM hindcasting applications are still scarce. Large carnivores and herbivores are archetypical examples of species whose current distributions are strongly determined by human disturbance ( Mladenoff , 1999 ; Kuemmerle , 2010 ). These species require large tracts of intact habitat, often conflict with land use, and are prone to poaching ( Breitenmoser, 1998 ; Gordon & Loison, 2009 ). As a result, large carnivores and herbivores are among the most challenging species to preserve in human‐dominated landscapes, and their populations have declined world‐wide. Safeguarding existing populations of large herbivores and carnivores and restoring their crucial ecological roles are thus important and well‐accepted conservation goals that may also benefit many other species ( Carroll , 2001 ; Vera , 2006 ; Gordon & Loison, 2009 ). Europe’s temperate zone has a long history of human settlement and land use ( Bramanti , 2009 ; Kaplan , 2009 ), and this has resulted in the extirpation of large carnivores and herbivores throughout the majority of their former ranges ( Breitenmoser, 1998 ; Vera , 2006 ). The European bison (or wisent, Bison bonasus L.) is Europe’s largest land mammal and the last surviving large grazer, as well as a prime example of a species that today only persists in a small, fragmented population. Bison disappeared from the wild in the early 20th century, but a systematic breeding and reintroduction programme prevented their extinction. Today, wild European bison occur in about 35 small, isolated herds in Central and Eastern Europe ( Pucek , 2004 ; Krasinska & Krasinski, 2007 ). As genetic diversity of the European bison population is low, the species continues to be at risk ( Tokarska , 2011 ), and the survival of wild bison will depend on establishing large metapopulations ( Perzanowski , 2004 ; Kuemmerle , 2011a ). However, identifying appropriate sites to establish large bison metapopulations is challenging because of two reasons. First, the habitat requirements of European bison are not fully understood. European bison are traditionally thought of as a species associated with closed temperate forests ( Heptner , 1961 ), but such ecosystems may simply have been their last refuges. Likewise, the species is generally considered a grazer, although browsing could contribute substantially to their diet ( Kowalczyk , 2011 ). Second, substantial uncertainty exists regarding both the species’ historic distribution ( Heptner , 1961 ; Pucek , 2004 ) and the effects of past environmental change and expanding human populations on the bison’s range. Recently, human population density and land use intensity have declined in many regions across Central and Eastern Europe in the wake of the collapse of socialism ( Ioffe & Nefedova, 2004 ; Baumann , 2011 ). This may offer a window of opportunity to implement a broad‐scale conservation strategy for safeguarding European bison and other large carnivores and herbivores ( Kuemmerle , 2011b ). A better understanding of habitat characteristics and historic distributions of European bison is thus urgently needed to identify priority sites for European bison conservation. Here, our goal was to reconstruct European bison range dynamics for the last 8000 years. To parameterize our SDM, we used a comprehensive set of Holocene European bison occurrences from before the species went extinct in the wild together with bioclimatic variables and vegetation reconstructions from a dynamic vegetation model as predictors. Including vegetation reconstructions addresses two limitations of SDM based solely on bioclimatic variables: First, bioclimatic variables tend to result in simplistic niche characterizations and may overpredict geographic distributions substantially ( Pearson & Dawson, 2003 ; Wiens , 2009 ). Second, climate‐based assessments of range shifts for species at higher trophic levels disregards that changes in vegetation communities, on which these species may depend, are slower than changes in climate. Dynamic vegetation models can capture such vegetation transitions ( Schurgers , 2009 ; Hickler , 2011 ). Using simulated vegetation in tandem with niche models can therefore improve range reconstructions for species at higher trophic levels, but to our knowledge, no prior study has done this. Specifically, we asked the following research questions: 1 What characterized the realized niche of European bison during the Holocene? 2 What was the geographic distribution of European bison during the last 8000 years and how did this distribution vary? 3 How have human population growth and farmland expansion affected the habitat of European bison? Methods Holocene records of European bison occurrence We used two types of historic European bison occurrences: (1) bone findings from archaeological sites and natural deposits, and (2) written records of bison occurrences. Locations of bison bones stemmed from a comprehensive archaeozoological database containing more than 7000 faunal assemblages from 31 countries across Europe since the last glaciation ( Benecke, 1999 ). This database contains 169 geolocated and dated findings of European bison remains ( Benecke, 2005 ), mainly from the mid‐ and late Holocene ( Fig. 1 ). Because the Benecke (1999) database is particularly rich in faunal assemblages from Southern, Western, Northern and Central Europe, but relatively sparse in Eastern Europe, we also included 45 occurrence locations from Heptner (1961) covering the European region of the former Soviet Union ( Fig. 1 ). These occurrence locations were derived from written records from various historical sources from the Middle Ages until the 20th century ( Heptner , 1961 ). Following Benecke (2005) , we assigned all occurrence locations to one of three time periods: (1) mid‐Holocene (6000 bc –1000 bc ), (2) Iron Age (1000 bc –600 ad ) and (3) Middle Ages and Modern Period (600 ad –1900 ad ). 1 European bison occurrences during the mid‐ and late Holocene based on bone findings ( Benecke, 2005 ) and written records ( Heptner , 1961 ). Only bison occurrences before the extirpation of the species in the wild in the early 20th century were used. The archaeozoological database ( Benecke, 1999 ) contains numerous faunal assemblages from Western and Southern Europe, but bison remains were absent from these regions during the mid‐ and late Holocene. (Projection: Albers equal‐area conic). Predictor variables Our study region covered all of Europe (15° W – 65° E and 75° N – 30° N). Climate variables for the dynamic vegetation model and the range reconstructions were derived from a transient simulation from 7000 bc to 2000 ad ( Schurgers , 2006 ) with a complex earth system model that includes atmosphere and ocean dynamics, as well as terrestrial and marine carbon cycling ( Mikolajewicz , 2007 ). The anomalies (averaged per 100‐year period) from this simulation compared to a 10,000‐year control simulation were superimposed on a detrended climate dataset with monthly average temperature, precipitation and cloud cover for the 20th century ( Mitchell & Jones, 2005 ). We then derived annual precipitation, potential evapotranspiration, growing degree days (temperature sum above a base temperature of 5°), as well as the minimum, mean and maximum temperature of the last year of each 100‐year period. Winter severity is a critical factor affecting ungulate survival, including that of European bison ( Krasinska & Krasinski, 2007 ). As a proxy variable for winter severity, we calculated the mean temperature of the coldest quarter (December–February). Vegetation reconstructions were simulated using an updated version of the dynamic vegetation model LPJ‐GUESS (v2.1, Smith , 2001 ). Vegetation was represented by 11 global plant functional types (PFTs) with different bioclimatic niches, growth forms, leaf phenology type (i.e. evergreen, summergreen or raingreen) and photosynthetic pathway (C3 or C4). The LPJ‐GUESS model successfully reproduces the main features of the global distribution of vegetation types, in particular in the study area. In this study, the model was run continuously for 10,000 years, using the first 1000 years to spin up the vegetation from bare ground. Model inputs consisted of the climate data described above, present‐day soil texture ( Sitch , 2003 ), and reconstructed and observed CO 2 concentrations ( Indermühle , 1999 ). As input variables for the SDM, we extracted the leaf area index (LAI; averaged per 100‐year period) for the PFTs that occurred in our study region. This yielded seven vegetation variables: LAI of boreal coniferous trees, LAI of temperate, broadleaved deciduous trees, LAI of temperate, shade‐intolerant deciduous trees, LAI of temperate evergreen trees, LAI of C3 grasses, total LAI and total woody LAI (see Fig. S1 in Supporting Information). All climate and vegetation variables were projected to the Albers equal‐area conic coordinate system. Reconstructing European bison range dynamics To characterize the realized niche of European bison during the Holocene, we used multitemporal niche calibration ( Nogues‐Bravo, 2009 ). This approach accounts for the fact that the environmental conditions a species utilizes may vary through time by matching each occurrence point with the environmental data from the time when a species occurred at a location. Using the resulting multitemporal dataset, we parameterized a single niche model that was then projected onto the environmental layers of each time period. Thus, all information about environmental conditions a species utilized throughout the full time span considered (6000 bc –2000 ad ) was used to describe its geographic distribution at a particular point in time ( Nogues‐Bravo, 2009 ). Calibrating niche models in this way results in more accurate and robust range reconstructions than reconstructions based on either past or present environmental predictors alone ( Broennimann & Guisan, 2008 ; Nogues‐Bravo, 2009 ). We extracted environmental conditions for the bison occurrences using the climate and vegetation layers closest to the mid‐point of our time periods (3000 bc , 0 ad , 1000 ad for the mid‐Holocene, Iron Age, and Middle Ages, respectively). We used maximum entropy modelling ( Phillips , 2006 ), to reconstruct the Holocene range of European bison. The maximum entropy approach assumes that the potential geographic distribution of a species is a probability distribution π over all locations (i.e. cells). To approximate π, a probability distribution is derived, considering constraints imposed by the environmental predictors at the occurrence locations ( Elith , 2011 ). The distribution with maximum entropy approximates π best because it is least constrained. Maximum entropy modelling requires only presence data, works well with small sample sizes and is among the highest performing SDM approaches ( Elith , 2006 ). A detailed mathematical description of maximum entropy modelling is provided in Phillips (2006) and Elith (2011) . To fit maximum entropy models, we used Maxent (v3.3.3). All model runs used default regularization, a maximum of 2500 iterations and 10,000 random background points ( Phillips & Dudik, 2008 ). Background points were randomly sampled from the environmental layers from all three time periods of our occurrence data. To prevent overfitting, we used only quadratic and hinge features. We built two groups of Maxent models: (1) models based on climate variables only and (2) models based on climate and vegetation variables. Furthermore, a few of our predictor variables were collinear ( r > 0.7) and we fitted alternative models in such cases, retaining the variable yielding higher goodness‐of‐fit. We assessed the goodness‐of‐fit of our models based on the receiver operating characteristics (ROC), which compares the true positive rate (i.e. sensitivity) versus the false positive rates (i.e. 1−specificity) across all possible thresholds, and by calculating the area under the curve (AUC), which provides a threshold‐independent measure of model performance ( Phillips , 2006 ). We used a fivefold cross‐validation strategy (i.e. five model runs based on 80% of the occurrence data, while retaining 20% for validation) and calculated mean AUC values and AUC standard errors for each model run. Final range maps were calculated as the average of the five replicate runs per millennia, and we used a logistic link function to yield a relative environmental suitability index (ESI) between zero and one ( Phillips & Dudik, 2008 ). We also used a jackknife procedure to assess variable importance, where we fitted models excluding a particular variable as well as single‐variable models and compared changes in AUC compared to the full model. Once Maxent models for both model groups were parameterized, we projected bison ranges for each millennia from 6000 bc to 2000 ad . Because European bison are a non‐migratory species ( Krasinska & Krasinski, 2007 ), our range maps depict the potential year‐round occurrence of European bison. To summarize the range maps, we extracted the area covered by ESI values larger than the minimum and the 10‐percent presence ESI values at our occurrence (i.e. training) locations for each millennium, and we calculated the median‐filtered average ESI map across all millennia. Reconstructing habitat fragmentation We used the HYDE 3.1 database that provides human population and land use grids for the last 12,000 years at a spatial resolution of 5 arc‐min ( Klein Goldewijk , 2010a,b ). We calculated the relative share of cropland and pasture for each grid cell and used human population density for all years for which ESI maps had been calculated (every millennium since 6000 bc ) plus the years 1800 ad and 1900 ad . All maps were transformed to the Albers equal‐area conic projection and resampled to 10‐km resolution. To assess how population and land use expansion affected the bison’s range, we crosstabulated the cropland, pasture and population density maps with the ESI maps. We used 5%‐wide bins for the cropland and grassland maps, logarithmic bins for the population density map (using the breaks 1, 10, 100, 1000, 10,000 and 100,000), and the minimum and 10% presence thresholds for the ESI maps. For 1800 ad and 1900 ad , we used the 2000 ad ESI map. To assess how human population and land use expansion fragmented the bison’s range, we masked all areas above the cropland, pasture and human population thresholds that were assumed to exclude bison populations. We compared eight scenarios: two levels of ESI thresholds (minimum and 10% presence) and four levels of land use and population thresholds: (1) cropland or pasture density > 5%/population density > 5 persons km −2 , (2) > 10%/> 10 persons km −2 , (3) > 15%/> 20 persons km −2 , and (4) > 20%/> 40 persons km −2 . Next, we derived the remaining habitat area, the number of habitat patches, mean patch size, the largest patch index (i.e. percentage of the landscape comprised by the largest habitat patch), and the mean proximity index, which measures isolation of patches ( McGarigal , 2002 ). As a neighbourhood distance for the proximity index, we used 300 km (i.e. roughly the maximum recorded dispersal distance of European bison, Krasinska & Krasinski, 2007 ; Kuemmerle , 2011a ). Results European bison inhabited a broad range of environmental conditions during the last 8000 years ( Fig. 2 ). Areas inhabited by bison had a range of average annual temperature between 0 and 13.2 °C, and bison tolerated average winter (i.e. December–February) temperatures as low as −14.9 °C. Bison also tolerated a wide range of annual precipitation, but avoided dryer regions (e.g. steppe regions). Bison occurrences were characterized by high vegetation cover and a mixture of woody vegetation and grasses. We did not find a clear preference for broadleaved versus coniferous forests, but bison occurred mainly in forests consisting of shade‐tolerant species ( Fig. 2 ). 2 Range of climate and vegetation conditions covered by bison occurrence locations during the mid‐ and late Holocene. Climate and vegetation information for the individual bison occurrence points (dated bone findings or written records) were extracted from the time periods when bison occurred at these locations and pooled to allow for a comprehensive description of the bison’s realized niche for the time period 6000 bp –2000 ad . Our final Maxent model induced three climate variables (average annual temperature, average winter temperature, and annual precipitation) and six vegetation variables (LAIs of boreal coniferous trees, temperate broadleaved trees, temperate shade‐intolerant broadleaved trees, temperate evergreen trees, C3 grasses and total LAI). The goodness‐of‐fit of this model was high, with an average cross‐validated AUC value of 0.920 and a standard error of 0.012. Based on our jackknife analysis, LAI of temperate broadleaved and LAI of shade‐intolerant broadleaved trees were the most important variables in our model, followed by mean winter temperature (which also contained the most information not present in any other variable). Variable response curves confirmed that environmental suitability for European bison was high in areas of high total vegetation cover that were dominated by temperate broadleaved forests but also contained grasses. Our models improved when including vegetation predictors compared to using climate predictors only. Our best model based on climate predictors alone included average annual temperature, average winter temperature and annual precipitation (AUC = 0.901). Distributional patterns differed substantially between the two model types, generally resulting in fewer and more clustered areas of high ESI values when including also vegetation predictors. For example, extensive areas in southern Scandinavia and Turkey appeared climatically suitable for European bison, but had low ESI values because of their vegetation composition (see Fig. S2). Projecting our calibrated niche model to all nine millennia and calculating average ESI values revealed that the heartland of European bison during the mid‐ and late Holocene was in Central and Eastern Europe with highest suitability values in contemporary Ukraine, Poland, Belarus, the Baltics, Romania and European Russia. ESI values were also moderately high in southern Scandinavia, the Balkans, parts of Turkey, and the Caucasus. In contrast, Western Europe contained few suitable areas ( Fig. 3 ). 3 Distribution of European bison during the mid‐ and late Holocene [average environmental suitability index (ESI) from 6000 bc –2000 ad ]. ESI values were derived by projecting a time‐calibrated niche model based on Holocene bison occurrences and a maximum entropy model to the environmental layers (climate and vegetation reconstructions) for each millennium. (Projection: Albers equal‐area conic). The geographic distribution of suitable areas was relatively stable during the last 8000 years ( Fig. 4 ). Variability in time occurred mainly in European Russia and southern Scandinavia, with a more northern distribution before 1000 bc . Today’s distributional patters are most similar to those from the mid‐Holocene (6000 bc and 5000 bc ). The distribution of high ESI values was not contiguous, but limited to isolated patches occurring, for example, in southern France and the Caucasus ( Fig. 4 ). Some of these patches, for instance the Caucasus, were not always isolated though, and the Caucasus did harbour a subspecies of European bison until the early 20th century ( Krasinska & Krasinski, 2007 ). The area covered by ESI values above the minimum and 10% presence thresholds did also not vary substantially (2.63 and 5.85 million km 2 on average, respectively, Fig. 5 ). 4 Dynamics in the distribution of environmentally areas suitable for European bison during the Holocene based on projecting a time‐calibrated niche model to the environmental layers of each time period. (Projection: Albers equal‐area conic). 5 Changes in the area characterized by environmental suitability index (ESI) values larger than the minimum value and the 10%‐percentile value observed at any of the locations where European bison were present during the mid‐ and late Holocene. Comparing our distributional maps with maps of cropland, pasture, and human population density suggests that the area suitable and available for European bison diminished steadily during the last 8000 years. For example, of the entire area with ESI values above the minimum presence threshold, almost 60% contained cropland ( Fig. 6a ) and more than 70% had a population density exceeding 10 persons km −2 ( Fig. 6e ) by 0 ad . Cropland expansion affected the bison’s range more than pasture expansion, although most areas with high ESI values contained pasture land by 1000 ad ( Fig. 6c,d ). After 1800 ad , areas with high ESI and low cropland, pasture and human population density decreased dramatically. For instance, of all the areas with ESI values above the 10% presence threshold, only 1.5% were without farmland in 2000 ad ( Fig. 6b ) and only 0.22% had human population densities < 10 persons km −2 ( Fig. 6f ). 6 Distribution of environmentally suitable areas for European bison among areas with varying cropland (a and b), pasture (c and d) and population densities (e and f). Left column: 100% = all areas with environmental suitability index (ESI) values larger than the minimum presence threshold; right column: 100% = all areas with ESI larger than the 10%‐presence threshold. Land use and human population expansion gradually fragmented the area suitable for European bison until about 1000 ad , after which fragmentation increased sharply ( Fig. 7 ). For example, the number of habitat patches increased steadily between 6000 bc and 1000 ad ( Fig. 7a ), while mean patch area ( Fig. 7b ) and largest patch index ( Fig. 7c ) decreased, indicating fragmentation. Range fragmentation reached a threshold in 1000 ad , whereafter all fragmentation parameters changed drastically. Our sensitivity analyses revealed that this general pattern did not depend on a particular farmland or human population density threshold (see Fig. S3). Bison habitat was more fragmented and less widespread when applying the 10% presence compared to the minimum presence ESI threshold ( Fig. 7 ). 7 Number of habitat patches, mean patch area, largest patch index and mean proximity index of the remaining European bison habitat after masking areas with high cropland, pasture or human population density. We masked areas with cropland or pasture density > 10% or human population density > 10 persons km −2 and used two ESI thresholds to derive habitat for each time layer: minimum presence (dashed line) and 10% presence (straight line) ESI values. Discussion Using a time‐calibrated niche model in tandem with a dynamic vegetation model, we found that European bison habitat preferences may be broader than previously assumed. During the mid‐ and late Holocene, the heartland of European bison was in Central and Eastern Europe, and our models suggested that bison had a more eastern and northern distribution than assumed by expert‐based assessments. Whereas climate and vegetation transformations did not affect the distribution of European bison substantially during the last 8000 years, human population growth and the expansion of farming drastically diminished and fragmented European bison habitat, likely reaching a tipping point during the last 1000 years. Our findings have important conservation implications, suggesting that efforts to preserve bison should move beyond the paradigm that bison are a species of broadleaved forests only and focus on the northern and eastern edge of its current distribution to establish large bison metapopulations. More generally, our study showed that hindcasting approaches can help to better characterize species niches and to reconstruct range dynamics, which is crucial to design effective conservation strategies that cover the entire ranges of endangered species. Traditionally, bison were considered to thrive only in interior temperate broadleaved forests, but this may simply be because such habitats were the species’ last refuges ( Pucek , 2004 ; Krasinska & Krasinski, 2007 ; Kerley , 2011 ). Our analyses showed that bison may utilize a wide range of forest types and occupy temperate broadleaved and southern boreal mixed forests alike. These findings are in line with both field studies ( Perzanowski , 2008 ), fine‐scale habitat assessments ( Kuemmerle , 2010, 2011b ) and expert‐based assessments ( Benecke, 2005 ). Likewise, the diet preference of European bison is not fully understood. Previous studies found browsing to play a minor role ( Gebczynska , 1991 ; Krasinska & Krasinski, 2007 ), and physiological studies suggest that bison are mainly grazers ( Mendoza & Palmqvist, 2008 ). Yet, our models showed that bison selected for areas containing both forest and grass cover, but were most widespread during the mid‐ and late Holocene in closed forests (as opposed to more open forests, Fig. 2 ). This suggests that browsing may have accounted for a substantial share of bison’s diet, confirming a recent field study finding a browsing share of up to 65% if supplementary feeding is excluded ( Kowalczyk , 2011 ). While it was interesting to see that we did not find support for bison inhabiting grass‐dominated regions (e.g. steppes in contemporary southern Russia and parts of Ukraine) during the mid‐ and late Holocene, we would like to point out that we assessed only the bison’s realized niche. The species’ fundamental niche could be wider, possibly including grassland conditions, especially when considering that European bison evolved from the steppe bison ( Bison priscus ) ( Pucek , 2004 ; Benecke, 2005 ). Likewise, although vegetation dynamics during the mid‐ and late Holocene did not affect bison distribution substantially, the transformation from grasslands to forest‐dominated ecosystems during the early Holocene likely affected bison distribution ( Pucek , 2004 ; Benecke, 2005 ). The Holocene distribution of European bison that our model predicted differed substantially from previous, expert‐based assessments of the bison’s range ( Heptner , 1961 ; Pucek , 2004 ; Benecke, 2005 ; Sipko, 2009 ; Tokarska , 2011 ). We found the heartland of European bison to be in Central and Eastern Europe (i.e. consistent with prior assessments), but the range of European bison varied strongly on its eastern and northern edge during the last 8000 years. These are also the regions where expert‐based range assessments disagree strongest. Overall, our analyses suggested that suitable regions for European bison stretched further north and east than previously thought, at least during parts of the Holocene. Bison’s range in these regions was limited by winter severity, the main factor determining ungulates survival in general and the most important climate variable in our models, which supports the ecological realism of our model. We also note that there is evidence that bison may have indeed inhabited these more northern regions ( Sipko, 2009 ). European bison distribution during the Holocene did not extend substantially into Western Europe (e.g. France or northern Spain). While this contrasts with most prior range assessments (but see Tokarska , 2011 ), and although the species was present there during the Late Stone Age (e.g. as documented in the cave paintings of Altamira and Lascaux) and possibly also the early Holocene, no bison remains from the mid‐ and late Holocene have been found there despite a large number of archaeozoological assemblages ( Vigne , 2003 ; Benecke, 2005 ). Some of these regions were grass‐dominated (e.g. Spain), had a high share of evergreen broadleaved trees (e.g. southern France), or forests with very high LAIs (and low grass cover), all of which differed substantially from the environmental conditions at our occurrence locations ( Fig. 2 ). Overall, this suggests that European bison were rare in Western and South‐western Europe already by the mid‐Holocene. We caution, however, that the reasons for this remain somewhat unclear, as land use and human disturbance may also have pushed bison out of otherwise suitable areas on the south‐western and western edge of their range (see below). Potential causes of the extirpation of European bison from large parts of their range are changing environmental conditions, habitat fragmentation and overexploitation. Environmental change has triggered range contractions and population collapses of large mammals in Northern Eurasia during the Holocene ( Schmölcke & Zachos, 2005 ; Nogues‐Bravo , 2008 ), but our results confirm that human factors were the main cause of the decline of European bison. Areas suitable for European bison were relatively stable both in geographic extent and location. On the other hand, the expansion of settlements and farmland during the Holocene severely reduced and fragmented European bison habitat, particularly in Southern and Central Europe. The expansion of human populations and farming may also explain why some suitable regions within the bison’s range were apparently not occupied during the mid‐ and late Holocene, such as the Balkans, Bulgaria, southern Romania or southern France. Settlements and farming were already widespread in these regions during the mid‐Holocene ( Bramanti , 2009 ), and forest cover was likely drastically reduced by 1000 bc ( Kaplan , 2009 ). On the other hand, later farmland expansion and lower human density may be the reason why bison could hold out longest in the north of Central and Eastern Europe ( Heptner , 1961 ), a region that maintained high forest cover until the 19th century ( Kaplan , 2009 ). Habitat loss and fragmentation did not occur uniformly in time, though. Both the extent and the connectivity of habitat plummeted during the last 1000 years, particularly since 1800 ad ( Fig. 7 ), when habitat availability was well below thresholds that may result in extirpation ( Fahrig, 2002 ). This suggests that the European bison population as a whole may have crossed a tipping point towards population collapse. Habitat fragmentation was also exacerbated by increased hunting, and overexploitation, which undoubtedly contributed to the decline and eventual extirpation of European bison. Our models yielded high goodness‐of‐fit and plausible distributional patterns that are consistent with independent assessments of the range contraction and population decline of European bison during the Holocene. A few sources of uncertainty remain. First, while we used a large set of occurrence data from both bone findings and written records, we cannot fully rule out bias in our occurrence data (e.g. because bone preservation may differ among soil types, sampling effort may be unequally distributed, B. bonasus remains cannot be clearly identified or bison did not occupy all environmentally suitable regions, Schmölcke & Zachos, 2005 ). While faunal assemblages were also widespread in regions where no bison records were found (e.g. France, Spain), bone findings were often associated with human settlements, adding uncertainty in areas that were settled later ( Bramanti , 2009 ; Kaplan , 2009 ). However, sampling bias does not affect our range reconstructions as long as environmental analogues (in either time or space) to undersampled regions exist in our occurrence dataset, representing a major advantage of our time‐calibrated niche model over mono‐temporal SDM or traditional range reconstructions based on visual interpretation of historic occurrence locations. Second, we used only a single climate and vegetation reconstruction. However, climate variations over this period are primarily driven by well understood orbital changes, and the main determinants of vegetation distribution and structure in Europe are well known and covered by LPJ‐GUESS ( Hickler , 2011 ). We also note that our vegetation constructions compare favourably to pollen reconstructions ( Prentice , 1998 ; Tarasov , 1998 ) and that we are not aware of any alternative climate reconstructions from Global Climate Models or transient vegetation reconstructions for the time period we assessed. Third, while our cross‐validation strategy allows for robust comparisons among alternative models, it may be optimistic in assessing overall model fit ( Araújo , 2005 ). However, cross‐validation remains the only feasible approach to assess the goodness‐of‐fit of our models as an independent set of presence/absence data cannot be gathered retrospectively ( Araújo , 2005 ). Last, while our sensitivity analyses showed that our main conclusions remain unaffected by the particular choice of land use and human population thresholds, the HYDE database has some uncertainties, especially regarding pastures and for older time layers ( Klein Goldewijk , 2010a,b ). Comparisons with other land use reconstructions suggest relatively consistent farmland expansion and deforestation patterns ( Olofsson & Hickler, 2008 ; Kaplan , 2009 ; Gaillard , 2011 ). Yet, the HYDE scenario likely represents a conservative estimate of the magnitude of land use change ( Kaplan , 2009 ), meaning that range fragmentation could have occurred even earlier and more rapidly than indicated by our results. Our range results have a number of important implications for the conservation of European bison and other large herbivores and carnivores. First, our results showed that the habitat preferences of European bison during the Holocene were broader than previously thought, with bison thriving in semi‐open areas and in broadleaved, mixed and coniferous forests alike ( Fig. 2 ). This provides further support for views that bison conservation should also focus on areas outside deciduous broadleaved forests and that the area of potential suitable European bison habitat may be large ( Pucek , 2004 ; Kerley , 2011 ; Kuemmerle , 2011b ). Second, our hindcasting approach suggests that European bison had a larger range in Europe’s north and east than previously assumed, yet that bison were absent or only present in very low numbers in Western Europe already during the mid‐ and late Holocene. As preserving European bison in the wild will depend on establishing large metapopulation ( Pucek , 2004 ), efforts to establish such metapopulations should focus especially on central European Russia and the Baltic states. Rural populations and land use there have declined drastically there since 1991 ( Ioffe , 2004 ; Baumann , 2011 ), which could benefit the establishment of a large bison population. Moreover, whereas unmanaged grasslands are sparse in Western Europe, abandoned and fallow farmland is now widespread in Europe’s East, allowing bison to utilize the broad range of habitat types they have used in the past. Third, future reintroductions should also focus on the northern edge of the current European bison distribution. Our range reconstructions showed that bison extended substantially more northwards in the past, especially during warmer periods ( Fig. 4 ). Considering that climate change will shift northern European regions to more temperate conditions during the 21st century ( Hickler , 2011 ), northern European Russia may offer suitable candidate sites for establishing a large bison metapopulation ( Sipko, 2009 ; Kuemmerle , 2011b ). Fourth, our study showed that drastic habitat fragmentation likely contributed substantially to the decline of the European bison. As linking isolated bison populations is a conservation priority identified in the European bison species action plan ( Pucek , 2004 ), this provides another argument for focusing European bison conservation efforts on Belarus, European Russia and the Baltics, where habitat fragmentation is still much lower than in Western and Central Europe ( Angelstam , 1997 ). Fifth, conservation managers should interpret low ESI values in parts of Western Europe carefully, as there is a possibility that bison were pushed out of suitable habitat there prior to our assessment period. However, even if these regions contained suitable habitat in the past, we believe that the prospects for the large bison metapopulations needed to ensure the persistence of the species in the wild are low because of conflicts with land use ( Kuemmerle ., 2011b ). Finally, our assessments confirmed that during the last 8000 years, bison in the Caucasus were separated several times from the rest of the European bison range, suggesting that the two bison lines may represent different subspecies and that conservation managers should continue to keep these lines separated. Reconstructing range dynamics of European bison suggested that land use change was the major force contributing to the catastrophic decline of European bison, whereas environmental change did not affect the distribution of bison substantially during the last 8000 years. Most future range assessments currently focus exclusively on possible future climate effects. Our study thus emphasizes the need to consider land use change when assessing species’ ranges and, more generally, the future of biodiversity ( Brook , 2008 ; Leadley , 2010 ). While conservation should not be about mimicking the past, hindcasting can improve our understanding of species’ niches and historic distributions. Such baselines are urgently needed for range‐wide conservation planning and to help conservation managers to design resilient strategies to preserve species in an increasingly human‐dominated world. Acknowledgements We thank L. Baskin, W. Cramer, C. Fløjgaard, P. Daskiewicz, K. Perzanowski, T. Samojlik, T. Sipko, J.‐D. Vigne and D.M. Waller for valuable discussions, and N. Benecke and V.G. Heptner for compiling the comprehensive datasets of historic European bison occurrences. Associate editor N. Roura‐Pasqual, J.P.G.M. Cromsigt and an anonymous reviewer are thanked for constructive comments on prior manuscript versions. We gratefully acknowledge support by the Alexander von Humboldt Foundation, the EU project ECOCHANGE (FP6‐036866), the NASA Biodiversity and NASA Land Cover and Land Use Change programmes, and the LOEWE initiative for scientific and economic excellence of the German federal state of Hesse. This study is a contribution to the Swedish Strategic Research area Biodiversity and Ecosystem services in a Changing Climate (BECC).

Journal

Diversity and DistributionsWiley

Published: Jan 1, 2012

There are no references for this article.