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

Learn More →

Modeling bird species distribution change in fire prone Mediterranean landscapes: incorporating species dispersal and landscape dynamics

Modeling bird species distribution change in fire prone Mediterranean landscapes: incorporating... Understanding species distribution dynamics is a central topic in ecology ( Brown and Lomolino 1998 ), and is of particular importance in the context of current rapid environmental change ( Thomas and Lennon 1999 , Thomas et al. 2004 ). The need to investigate and, if possible, anticipate the ecological impacts of current global environmental change, has led to the development of several modeling approaches for predicting future biodiversity spatial patterns ( Huntley et al. 2007 , Jetz et al. 2007 , La Sorte et al. 2009 ). The predictive modeling of species distributions based on niche theory has received a great deal of attention, and has significantly contributed to the quantification of the impacts of global change, such as climate change or the potential impact of invasive species on ecosystems ( Guisan and Thuiller 2007 , Thuiller et al. 2008 , Elith and Leathwick 2009 ). Species distribution models (SDMs) are founded, however, on the assumption that species distributions are in equilibrium with environmental conditions and that these species environment relationships can be used to estimate species responses to environmental changes ( Elith and Leathwick 2009 ). Relaxing this equilibrium assumption has been identified lately as key challenge in the field of species distribution modeling ( Dormann 2007 , Elith and Leathwick 2009 ). One way to begin addressing this challenge is to separate potential habitat estimation from the realized species distribution ( Soberon 2007 , Soberon and Nakamura 2009 ). Under non‐equilibrium conditions, habitat availability (areas in which environmental conditions allow species presence) and occupancy or abundance patterns differ ( Brotons et al. 2004 ). Recent studies have adopted hybrid SDM approaches combining correlative, phenomenological models with more mechanistic models, thereby allowing for the integration of critical ecological processes into SDMs ( Gallien et al. 2010 ). These new approaches purport to be more ecologically sound and achieve a better fit to the observed dynamics in species distributions ( Keith et al. 2008 , Elith and Leathwick 2009 , Willis et al. 2009 ). Such hybrid models explicitly incorporate ecological processes other than habitat selection (e.g. movement or competition), to better reflect actual species distribution patterns. Hybrid models allow us therefore to consider the dynamic processes of species distribution needed to model the effect of global change. As habitat connectivity and species dispersal capabilities are key factors determining a species use of newly appearing habitat and, hence, driving species responses to changes in habitat quality or availability it is important to include these processes in hybrid models ( Hanski and Ovaskainen 2000 ). However, dispersal has been explicitly incorporated in SDMs only in the last years ( Willis et al. 2009 ), and at mostly it has been used to trim model outputs using a rather static approach ( Peterson et al. 2001 , Allouche et al. 2008 ). These dynamical hybrid models will require however using datasets collected through time from validation purpose ( Keith et al. 2008 , Vallecillo et al. 2009 ). Here, we compare the ability of correlative and hybrid SDMs to quantify the distribution dynamics of an early succession bird, the ortolan bunting Emberiza hortulana , in a highly dynamic landscape driven by land abandonment and wildfire in Catalonia, north‐eastern Spain. We describe changes in this landscape using a landscape dynamics model that simulates fire impacts in the landscape and forest regeneration. Our hybrid model integrates the tracking of habitat suitability changes from a correlative model with a spatially‐explicit population model that incorporates dispersal, and assumes equilibrium in the initial state only. We investigate alternative scenarios of dispersal capability to generate predictions of distributions in the context of change in habitat suitability driven by fire disturbance. We finally validated correlative and hybrid model results using two independent data sets: one obtained from monitoring efforts deployed in recent burns; and the other obtained from a general‐purpose monitoring program, thus allowing models to be evaluated using a temporal perspective. The comparison of correlative and hybrid model performance allows us to test the degree to which species distribution dynamics (i.e. colonization patterns) are constrained by dispersal and/or habitat availability and thus evaluate the expected gain in model performance achieved by moving from purely correlative to hybrid SDMs in the case of species inhabiting highly‐dynamic landscapes. Methods Bird data The ortolan bunting is an open habitat, early‐succession bird species typically found in montane and subalpine meadows, open shrublands, as well as farmland, with a strong preference for recently burnt areas ( Brotons et al. 2008 , Menz et al. 2009 ). However, the analyses of large scale atlas data suggest that the use of recent burns may be constrained by dispersal thus making this species a good candidate to investigate species distributional changes in non‐equilibrium contexts ( Brotons et al. 2008 ). We compiled a range of large scale monitoring data relating to the distribution and spatial dynamics of the species in the region. First, comprehensive snapshot distribution data were obtained from the Catalan breeding bird atlas (CBBA) ( Estrada et al. 2004 ). The CBBA atlas sampled Catalonia (north‐eastern Spain) between 1999 and 2002 using a total of 3150 1 3 1 km plots in which presence/absence records were obtained from two one hour transects within each plot. Second, data on the species spatial dynamics were obtained from two bird monitoring programs from 2006 to 2009 in either unburnt (with relatively static habitat changes) or burnt (with more dynamic habitat change) sites. Information from unburnt sites was obtained from the Common Bird Monitoring Survey of Catalonia (CBS) ( Brotons et al. 2007 ). This general‐purpose monitoring program was initiated in 2002 and aims to track bird communities annually using 3 km transects divided in 500 m sections in which the number of observation of each species is noted at different distance bands. CBS transects are distributed within specific 10 3 10 km squares selected by volunteer observers throughout Catalonia. We compiled 3178 presence–absence records from CBS transects in 2006‐2009. Information from burnt sites was obtained from a similar monitoring method used to investigate bird species responses to recent wildfires in Catalonia (DINDIS) ( Zozaya et al. 2010 ). In this survey, started in 2006, all fires that occurred in Catalonia after year 2000 that are >50 ha are monitored yearly using a number of transects proportional to fire area. We compiled 1914 presence–absence records from DINDIS transects conducted between 2006 and 2009 on the 56 fires that occurred between 2000 and 2007. Environmental data and landscape dynamics Catalonia is a region with Mediterranean climate that is remarkably heterogeneous due to sharp climatic and geological gradients, including a range of landscapes from alpine habitats to coastal marshes, and from evergreen forests to steppes and agricultural mosaics. To describe habitat diversity relevant to our bird study model, we used geographical information system tools, to gather and combine information from available forestry, land cover and fire data ( Ministerio de Medio Ambiente 2006 ) to identify forested areas, their dominant species, shrubland and non‐forested open areas. The resulting forest map is composed of 50 3 50 m cells belonging to different forest (i.e. dominant species) and land use classes. Land abandonment and large wildfires appear to have driven regional landscape dynamics in recent years ( Díaz‐Delgado et al. 2004 , Gil‐Tena et al. 2009 ). We simulated temporal changes in landscape associated with these two processes using a spatially‐explicit landscape dynamic model implemented in the SELES (spatially explicit landscape event simulator) modeling platform ( Fall and Fall 2001 ). The landscape model includes a fire sub‐model and a shrub succession sub‐model. The fire sub‐model uses observed burnt areas between 2000 and 2009 directly and post‐fire vegetation changes were only modeled within these areas. Fire information was obtained from the historical sequence of fires available in Catalonia (Catalan Government, < http://mediambient.gencat.net/cat/el_medi/incendis/ >). Forest maturation was tracked using the years since the last fire. Post‐fire changes in forest composition (i.e. changes in dominant tree species or changes to shrub) and the probabilities for those transitions were obtained using bibliographic sources ( Rodrigo et al. 2004 ). We induced spatial autocorrelation in forest composition changes by coupling transitions of some cells with those of their neighborhood. Succession mimicking land abandonment impacts was only modeled for shrub (i.e. forests were assumed to be in a stable state). Shrub had a yearly probability of changing to forest after 20 yr. In case that a transition from shrub to forest occurred, the resulting forest type was chosen by using probabilities obtained from the distribution of mature forest types in the surroundings cells (150 m radius). The landscape simulation model was run 10 times for the period between years 2000 and 2009. Using the sequence of forest maps, produced by the landscape model, and additional topographic data, we then defined a set of landscape variables from which bird habitat suitability models could be derived (Supplementary material Appendix 1, Table A1). As the landscape model operates on a 50 3 50 m resolution grid, we averaged cell values for these landscape variables within a 1 3 1 km resolution grid to match the bird atlas and monitoring data. Correlative species models We developed three correlative SDMs of the distribution of ortolan buntings in the period 2000–2009 on a 1 3 1 km resolution grid ( Fig. 1 ). These models were based on different assumptions of species responses to environmental variables as follows. 1 Conceptual representations of the species distribution models used. The three correlative species distribution models; initial probability of occupancy (IPO); habitat suitability model (HS) and the combined habitat suitability model (CHS) were based statistical representations of habitat suitability, whereas the hybrid species distribution model included population dynamics and dispersal allowing habitat suitability to be tracked independently of occupancy. 1) Initial probability of occupancy (IPO) This model hypothesizes a static distribution throughout the studied period based on initial estimates of the species probability of occupancy based on the probability of occupancy for ortolan buntings in year 2000 ( Fig. 2 ), taken from ( Brotons et al. 2007 ). This correlative SDM included both environmental information as well as the spatial structure of bird distribution derived from the CBBA data (i.e. using autocovariates or other spatial predictors in the model) ( Augustin et al. 1996 , Lichstein et al. 2002 ). We assume that the IPO provides a good approximation for habitat suitability at year 2000. 2 Initial probability of occurrence (IPO) for ortolan bunting in Catalonia, taken from Brotons et al. (2007). 2) Pure environmental habitat suitability (HS) Given the hypothesis that habitat selection in the ortolan bunting is at least partially driven by changes in open habitats related to fire ( Menz et al. 2009 ), we used land use and forest variables describing burn areas and open‐habitat availability as explanatory factors (Supplementary material Appendix 1, Table A1) in a generalized linear model with binomial response. Model selection was based on AIC. To avoid strong departures from equilibrium in the training data, we calibrated the HS model using only CBBA data sampled within 10 km of a recorded presence. The highest ranked model included variables representing the percentage of young forest (<20 yr), and the percentage of unburnt and burnt shrubland in the 1 3 1 km cell (Supplementary material Appendix 1, Table A2). As shrubland and forest variables are updated by the landscape dynamics model, the HS derived an expected probability of species presence following changes in the landscape resulting from fire disturbance and succession processes. 3) Combined habitat suitability (CHS) In this third model, the IPO and HS models were combined to represent the effect of landscape changes while preserving information regarding the initial probability of occupancy. That is, we used IPO as an approximation to initial habitat suitability (CHS 0 =IPO) and used changes in HS to indicate increases or decreases of local habitat suitability. This leads to the following expression for such combined habitat suitability (CHS) for each cell at time t >0: CHS t = CHS t−1 + (HS t − HS t−1 ) where HS t is the value of the pure environmental habitat suitability model at time t . To arithmetically operate with the values of the IPO and HS, both models were expressed as probabilities. In the event that IPO=HS 0 , the environmental model would account for the species spatial structure and the equation above would simplify to CHS t =HS t . Hybrid species distribution modeling Hybridization of phenomenological with mechanistic models ( Gallien et al. 2010 ) facilitates the estimation of effective occupancy, defined as areas in which the species is likely to be present given both habitat suitability and ecological constraints on distribution limiting occupancy (e.g. dispersal limitation, competition) or favoring use of non‐favorable habitat (e.g. spillover effects). Here, we hybridized the combined habitat suitability (CHS) model with a spatially explicit population model. Assuming initial equilibrium, the hybrid model modifies an initial occupancy state at a given site (IPO) by tracking local changes in HS and the occurrence of the species in the neighborhood ( Akçakaya et al. 2005 ). Decreases in habitat quality will result in a progressive decrease in the individuals’ carrying capacity at that location and eventually resulting in local extinction. Increases in habitat quality in unoccupied sites (i.e. disturbance effects) lead to increases in the carrying capacity at that location, but these areas can only be colonized from neighboring occupied cells (i.e. population sources). Dispersal is, therefore, implemented as a key ecological constraint and leads to a possible delay in the occupancy of suitable habitat. When dispersal constraints are strong and habitat changes are fast, non‐equilibrium situations can arise rapidly. Although other mechanistic processes governing distribution, such as species interactions ( Willis et al. 2009 ), could be included, previous studies identify dispersal as the main constraint on ortolan bunting distribution ( Brotons et al. 2008 ). Like the spatially explicit fire‐landscape model, the spatially explicit population model was implemented in the SELES modeling platform ( Fall and Fall 2001 ). The population model is an age‐structured stochastic population model that tracks occupancy, abundance and carrying capacity for ortolan bunting males within each cell. The population model requires several inputs data and parameters which are described below. Initial occupancy and habitat suitability Initial occupancy and habitat suitability patterns are stochastically determined using the probabilities indicated in the correlative models. For any given cell i , a random uniform number between 0 and 1 is drawn and this number is compared to the IPO i value to determine initial cell occupancy. The same random number is compared to CHS i t to determine cell suitability at any given time t . Using the same random number ensures that: 1) there are not occupied cells in unsuitable habitats at time 0 (because CHS 0 =IPO); 2) cells are only allowed to change their suitability when changes in CHS occur (i.e. if CHS value does not change from t to t+ 1, the suitability state at t is the same as the state at t+ 1). Initial abundance and changes in carrying capacity Within the initially occupied cells the number of males in the cell was determined from a Poisson distribution with a mean based on a step function relating the initial probability of occurrence with density data from the breeding bird survey of Catalonia (CBS) from the 20022004 period ( Table 1 ; Brotons et al. 2007 ). The initial number male juveniles is then determined using a constant adult:juvenile ratio. The carrying capacity ( K ) at any time t was set using the same step function ( Table 1 ), but with CHS t as input instead of IPO. Whenever a previously suitable cell becomes unsuitable then K is set to zero, and if the cell was occupied then its population was forced to disperse. 1 Parameterization of the spatially explicit population model (SEPM) for ortolan bunting. Parameter Value (s) A) Initialization (from Brotons et al. 2007 ) (0.01–0.02): 6 IPO: abundance relationships (0.02–0.07): 8 (occurrence probability): (number of individuals) (0.08–0.22): 15 (same for CHS‐K relationships) (0.23–1): 18 B) Demographic parameters (from Steifetten and Dale 2006 ) Initial adult/juvenile proportion 0.60 Breeding probability N(0.65, 0.0045) truncated to avoid values outside [0,1] Offspring mortality probability 0.146 Juvenile mortality probability 0.67 Adult mortality probability 0.375 Probability distribution of a clutch of size j ( j =0, 1, …, 4) [0, 0.02, 0.14, 0.43, 0.41] Male proportion at birth 0.50 C) Dispersal parameters (following Dale et al. 2005 , 2006 ) Proportion of dispersal for adults 0.50 Proportion of dispersal for juveniles 1 Parameter α for negative exponential dispersal kernel D1: 1; D2: 0.5; D3: 0.355; D4: 0.1 Local population dynamics We designed and parameterized mortality and reproduction for ortolan bunting according to the literature ( Steifetten and Dale 2006 ) and local data. The first local breeding occurs at the end of the first year (in autumn), followed by local mortality (applied in winter). Reproduction and mortality events occur once per year. The probability an individual reproduced varied stochastically following a truncated normal distribution and, if reproduction occurred, the clutch size was determined from a multinomial distribution ( Table 1 ). For each newborn, the number of male offspring increased by one if the sibling is male and does not die after birth. The number of male juveniles for a given year is set to the male offspring, and the surviving juveniles from the last breeding event become adults. Different mortality rates were applied to adults and juveniles, which determined the number of individuals surviving the winter in each occupied cell. Note that this number may be bigger than the local carrying capacity. Dispersal Individuals disperse into neighboring cells with adequate habitat according to a dispersal kernel function ( Chapman et al. 2007 ). In most scenarios (see below), we modeled the distribution of dispersal distances (i.e. the dispersal kernel) as a negative exponential function of geographical distance ( d ): e ad , where the inverse of parameter α can be interpreted as the mean dispersal distance. To speed up computations, we truncated all kernel values below 0.001 to 0. The first dispersal event occurs one year after the simulation starts, and subsequent dispersal events are scheduled for every summer thereafter. According to the available information on the species dispersal capabilities ( Dale et al. 2005 , 2006 ), the proportion of individuals dispersing from each cell is 50% for adults and 100% for juveniles ( Table 1 ). However, if remaining adults still exceed the current carrying capacity, the additional individuals are forced to migrate. The location to where individuals can spread is limited to suitable cells that can be reached according to the dispersal kernel function and where the population is below the carrying capacity. We assume that habitat type does not affect dispersal ability. When there are no suitable locations available the disperser is assumed to die due to the unavailability of resources in the surroundings of the source cell. We included an additional parameter in our model to account for behavioral constraints in dispersal induced by conspecific attraction leading to biased dispersal to areas where the species is abundant ( Boulinier and Danchin 1997 , Keitt et al. 2002 ). We modeled conspecific attraction by multiplying the dispersal kernel of occupied cells by a constant factor (see below). Dispersal scenarios and model runs We ran the population model under eleven different dispersal scenarios: no dispersal (ND) – all dispersed individuals stay in the source cell, and those exceeding the carrying capacity of the cell are assumed to die; restricted dispersal (D1D5) – in these scenarios we tested three different values of the negative exponential parameter. In scenario D1 (α=1, 1/α= 1 km) dispersers can move up to 7 km (that is, with α=1 kernel values are lower than 0.001 for distances >7 km and we truncated those kernel values to 0), in D2 (α=0.5, 1/α=2 km) up to 13 km and in D3 (α=0.355, 1/α=2.81 km) up to 18 km. The distribution of dispersal distances in D3 reflects that observed in the ortolan bunting population in Norway ( Dale et al. 2005 ); long dispersal (D4D5) – in these two scenarios dispersers could move up to 50 km. In scenario in D4 we still used the negative exponential kernel (α=0.1, 1/α=10 km), while in D5 the distribution of dispersal distances was uniform between 0 and 50 km; treatments with conspecific attraction (D1CD5C) – we replicated the previous scenarios (D1D5) by adding an effect of conspecific attraction to dispersal. Specifically, we multiplied by five the value of the dispersal kernel of cells that were already occupied, whereas on suitable but unoccupied cells the dispersal kernel value was not modified. The five‐fold increase was chosen based on a sensitivity analysis that showed weak effects for small values. The spatially explicit population model was run for the period of 2000 to 2009 with 100 replicates for each combination of landscape simulation and dispersal scenario, and replicates were averaged to obtain the overall probability of occupancy. Model runs and validation Correlative and hybrid SDM predictions were validated using the two monitoring datasets. Observed occupancy in non‐burnt sites (3178 presence–absence records) obtained from the CBS monitoring program primarily corresponds to areas with small decreases or no change in habitat suitability. Observed occupancy data (1914 presence–absence records) obtained in burnt sites from the DINDIS project corresponds to areas affected by large wild fires ( Zozaya et al. 2010 ), where habitat quality is predicted to strongly increase and eventually lead to species colonization ( Menz et al. 2009 ). We used two different validation approaches to assess its predictive accuracy. In the first approach we calculated the agreement between the modeled probability of occurrence and the observed occurrence using area under the receiver operating characteristic curve (AUC) statistics ( Fielding and Bell 1997 ). In the second approach (Supplementary material Appendix 1, Table A3) we estimated the log‐likelihood of the model probabilities (in a binomial distribution) given observed occurrence data. We also estimated the coefficient of determination for the model (R 2 ) by comparing its log‐likelihood value with the log‐likelihood of a null model where the species prevalence is the only information available ( Nagelkerke 1991 ). Validation statistics were computed for two periods (20062007 and 20082009) taking within each period the maximum observed and modeled values among the two years. As time progresses and fire impacts accumulate in the landscape, we expect a progressive divergence in model predictive ability between correlative approaches and hybrid models including species dispersal constraints. Results Landscape, habitat suitability and occupancy dynamics Given that in our landscape habitat quality was generally decreasing (i.e. widespread land abandonment leading to a loss of open habitats), the total amount of suitable cells in the HS model, and hence in the CHS model, generally decreased ( Fig. 3 ). Only small increases were observed in years where the general decrease in availability was compensated for by the creation of new habitat due to fires. In the scenario in which a species could not disperse from the initially occupied areas (ND in Fig. 3 ), the total amount of occupied cells decreased following decreases in habitat suitability due to encroachment, forest succession and ageing. In the high dispersal scenario (D5 in Fig. 3 ), new habitats created by fire could be reached and, therefore, allowed the occupancy of new patches leading to progressive shift in distribution from areas initially occupied, but losing quality, to these post‐fire newly created habitats. In the scenarios in which dispersal occurred and individuals had the capability to move and find new habitat, occupancy patterns were between the patterns showed by the two extreme behaviors described above (D1 and D5 in Fig. 3 ). The higher the dispersal capability, the closer the match between real distribution patterns and habitat suitability changes will be. When dispersal is constrained, species distribution showed non‐equilibrium situations and increasingly differed from habitat suitability estimates (D1 and CHS in Fig. 3 ). Furthermore, our results indicated that conspecific attraction acted by effectively constraining dispersal and limiting occupancy to newly appearing habitat (e.g. compare D1 and D1C in Fig. 3 ). 3 Predicted habitat suitability and occupancy dynamics for the hybrid model under different dispersal scenarios. The thick solid line represents the number of sites with potential habitat in the region (i.e. sum of CHS values) in each year. The thin solid line represents the predicted number of occupied cells in each year under a non‐dispersal scenario (ND), which indicates losses of occupancy associated with decreases in habitat quality. The remaining dashed lines represent occupancy dynamics under different scenarios of dispersal capability, with and without conspecific attraction (scenarios D1C, D3C and D5C and scenarios D1, D3 and D5 respectively). Dispersal capability allows the species to partially compensate for the loss of occupied sites by making use of new suitable habitat associated with fire impact. Model validation results Among all comparisons, the HS model performed worse than the CHS correlative model ( Fig. 4a and Supplementary material Appendix 1, A3), indicating that considering habitat suitability changes starting from an accurate initial species distribution was more effective than a pure environmental niche model. In cases in which species data were gathered from the monitoring of non‐burnt sites, we detected weak differences among dispersal scenarios of the hybrid model and a good overall performance of the IPO model ( Fig. 4a , Supplementary material Appendix 1, A3). Moreover, although all the models were satisfactory in terms of agreement with observed data, the models including distribution dynamics performed worse than the scenario assuming a static distribution (IPO). Thus, changes affecting habitats already occupied in the initial conditions were very limited and hence initial occupancy predictions fit monitoring observations adequately. 4 AUC evaluation results of the hybrid model under different dispersal scenarios, using either the (a) common bird survey data (CBS) or (b) monitoring data from burnt habitats (DINDIS). Correlative models (IPO, HS and CHS) are indicated using black symbols. Dispersal scenarios are grouped into five categories: no dispersal (ND), restricted dispersal (D1D3), restricted dispersal with conspecific attraction (D1CD3C), long dispersal (D4D5) and long dispersal with conspecific attraction ( D4C D5C ). Strong differences between dispersal scenarios occurred for DINDIS evaluation data (20062007: F 4,105 =176.6, p<1 e‐15; 20082009: F 4,105 =1141.4, p<1 e‐15), but not for the CBS evaluation data (20062007: F 4,105 =3.4, p=0.0195; 20082009: F 4,105 =2.0, p=0.1208). However, when we used monitoring data from burnt habitats, the scenario predictions showed more marked differences both in the likelihood of the observed data and the predictive power to capture the ortolan bunting presence and abundance patterns after fire impact ( Fig. 4b , Supplementary material Appendix 1, A3). The scenarios that were most successful at predicting which burned areas were colonized by buntings included restricted dispersal capabilities (scenarios D1D3 and D1CD3C); both the no dispersal (ND) and long dispersal (D4/D5 and D4C/D5C) scenarios were less successful at prediction ( Fig. 5 ). Furthermore, in accordance with our prediction, we observed that differences between models increased in time, with outputs from the hybrid SDM with restricted dispersal capabilities offering the best fit to observed data ( Fig. 4b ). Thus, ortolan buntings appeared to have the ability to colonize newly appearing suitable habitats, but these capabilities were limited by dispersal constraints, leading to non‐equilibrium situations between species distribution and habitat suitability. 5 Probability of occurrence on fire patches in 2009 according to the spatially explicit population model for dispersal scenario D1 and landscape (run no. 3). Bird symbols indicate those fire patches where ortolan bunting was observed in monitoring data on burnt sites (DINDIS project). Discussion We have shown for the ortolan bunting that dynamic species distribution modeling based on the explicit integration of habitat suitability and effective occupancy estimation may successfully capture species distribution dynamics even when species distributions are not in equilibrium with habitat conditions. This is achieved by including key ecological processes that appear to limit effective occupancy of suitable habitat (i.e. dispersal constraints). Our hybrid modeling approach facilitates the inclusion of such processes by using initial species distribution information and a careful estimation of potential colonizer sources, and therefore dispersal fluxes, into new appearing habitat. New habitat arises as a function of explicit habitat selection hypotheses coupled with a landscape dynamics model. Although hardly surprising, our models suggest that when changes in landscape suitability are abrupt, a hybrid approach, including dispersal limitations, is superior to purely correlative dynamic approaches that project habitat suitability under the new environmental conditions. The benefit of hybrid dynamic SDMs over purely correlative SDMs may be especially apparent in the cases in which validation data is used to assess dynamic distribution models ( Araújo et al. 2005 ). Regardless, models should be able to reproduce observed dynamics in areas with sharp changes in habitat suitability and not only overall distribution patterns dominated by areas with more stable occurrence patterns. In this context, we encourage the validation of dynamic distribution models using temporally replicated biological data ( Araújo et al. 2005 , Willis et al. 2009 ). Projections of distribution change that are not validated with temporally replicated data should be used with caution and interpreted as a general optimistic hypothesis of potential distribution change ( Gregory et al. 2009 ). Such projections have been previously applied to the investigations of bird responses to changes in fire regime using a metapopulation modeling approach linked to a landscape simulation model similar to ours ( Akçakaya et al. 2005 ). However promising, without independent temporal validation of the models it may become difficult to assess the relevance of the model output in terms of the validity of the assumptions specifically implemented into dynamic SDMs. To date, hybrid modeling approaches have been used primarily as tools to assess the impact of climate change on the viability of populations and species distribution shifts ( Wintle et al. 2005 , Brook et al. 2009 , Willis et al. 2009 ). Future applications of hybrid species distribution models should allow the study of interactions between land use changes, climate change and changes in the associated perturbation regimes. However, the capture of ultimate processes allowing the track of changes in habitat quality for a given species and their match with adequate landscape models is probably not an easy task. In our case, the capability of our habitat suitability model (HS) to assess changes in habitat quality was constrained by the capability of the landscape model to describe landscape changes. These limitations have probably led the hybrid model to perform worse on stable areas than the initial probability of occupancy (IPO), causing the IPO to rank better in the validation exercises. Therefore, careful integration of initial distribution conditions, key ecological processes and changes in habitat suitability across space are, therefore, critical challenges ahead. Our results suggest that these factors will all be important and vary in importance in line with the ecological context considered. For the ortolan bunting, suitable habitat arises as a result of the alternative processes of land abandonment, leading to decreases in habitat availability, and fire, leading to the creation of new habitat patches, making assessment of the overall impacts on species distribution complex ( Brotons et al. 2008 ). Early succession species may respond in very different ways to such complex shifting mosaics depending of their ability to find newly available habitats ( Brotons et al. 2005 ). Dispersal has been shown to be a key parameter in the distribution dynamics of a species in a context of changing habitat availability ( Hanski and Ovaskainen 2003 ). Our models indicate that the ortolan bunting has a limited capability to colonize new habitats and is, therefore, only partially able to compensate for habitat loss in undisturbed habitats by colonizing new burnt sites. We suggest that limited dispersal capability could be driven in part by conspecific attraction that may restrict effective dispersal distances ( Serrano and Tella 2003 ). Model agreement with observed data was highly sensitive to changes in dispersal distance (directly, or induced by conspecific attraction). These results suggests that previous models incorporating dispersal in a static way may be flawed as dispersal cannot be separated from the distribution context in which it takes place ( Peterson et al. 2001 , Allouche et al. 2008 ). Given the limited dispersal capability of the ortolan bunting, and the widespread loss of open habitat, the spatial pattern of the new habitat is predicted to be an important driver of species persistence. Such a sinking–emerging island model, where new habitat suddenly appears and gradually vanishes in time, can be useful for a number of other species in a context of rapid environmental change. In such cases, restricted distribution ranges coupled with strong dispersal constraints will likely lead to non‐equilibrium situations and high variability in the distribution outcome depending on the spatial pattern of habitat change. Positive spatial autocorrelation in fire distribution will likely favor metapopulation persistence in the long term, whereas new distant fires may make new colonization difficult thereby decreasing chances of metapopulation maintenance. The prediction of species distribution dynamics is contingent upon estimates of dispersal capabilities. Our results suggest that, if adequate data exist, hybrid species distribution models could be used to estimate model parameters accounting for the dispersal ability better fitting the observed species spatial dynamics ( Smolik et al. 2010 ). This reverse modeling approach would be of similar nature to the rationale used to parameterized metapopulation models ( Lele et al. 1998 , Étienne et al. 2004 ). Given the accumulation of temporal biological data and adequate information on species distribution dynamics (e.g. monitoring schemes for a number of taxa; Yoccoz et al. 2001 ), this kind of approach may allow the estimation of dispersal capabilities for a large number of species. The time and spatial scale of the training data must also be appropriate to the rate of landscape dynamics and allow the capture of enough distribution changes since otherwise the dispersal ability of the species may be underestimated. To some extent the success in ecology of metapopulation theory and species distribution modeling ( Guisan and Zimmermann 2000 , Hanski and Ovaskainen 2000 , Guisan and Thuiller 2007 ) is rooted in the use of widely available biological data to parameterize spatial or temporal species distribution dynamics. Hybrid modeling approaches linked to the use of large scale temporally replicated biological data sets have the potential to become a key tool in setting the scenarios of species responses to global change impacts (Thuiller et al. 2008, Elith and Leathwick 2009 ). Acknowledgements We greatly acknowledge the task of all the contributors to the projects run by the Inst. Català d’Ornitologia, whose careful work has been indispensable for the development of the present study. We also thank Svein Dale for information on the ortolan bunting dispersal data and the GRAF (Catalan fire fighters brigade) for providing information on fire perimeters. Both Catherine Graham (the editor) and Hawthorne Beyer provided valuable comments: we are in their debt. This work has received financial support from the projects CGL2005‐2000031/BOS, CGL2008‐05506‐C02‐01/BOS and Consolider Montes CSD2008‐00040 granted by the Spanish Ministry of Education and Science (MEC). The study is also a contribution to the FP7‐226852 project SCALES. Partial funding was also provided by the Catalan Government grant SGR2009‐531 and a Beatriu de Pinós contract to MC (2009 BP‐B 00342). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Ecography Wiley

Modeling bird species distribution change in fire prone Mediterranean landscapes: incorporating species dispersal and landscape dynamics

Download PDF
10 pages

Loading next page...
 
/lp/wiley/modeling-bird-species-distribution-change-in-fire-prone-mediterranean-AZ01QkrRIg

References (64)

Publisher
Wiley
Copyright
© 2011 The Authors
ISSN
0906-7590
eISSN
1600-0587
DOI
10.1111/j.1600-0587.2011.06878.x
Publisher site
See Article on Publisher Site

Abstract

Understanding species distribution dynamics is a central topic in ecology ( Brown and Lomolino 1998 ), and is of particular importance in the context of current rapid environmental change ( Thomas and Lennon 1999 , Thomas et al. 2004 ). The need to investigate and, if possible, anticipate the ecological impacts of current global environmental change, has led to the development of several modeling approaches for predicting future biodiversity spatial patterns ( Huntley et al. 2007 , Jetz et al. 2007 , La Sorte et al. 2009 ). The predictive modeling of species distributions based on niche theory has received a great deal of attention, and has significantly contributed to the quantification of the impacts of global change, such as climate change or the potential impact of invasive species on ecosystems ( Guisan and Thuiller 2007 , Thuiller et al. 2008 , Elith and Leathwick 2009 ). Species distribution models (SDMs) are founded, however, on the assumption that species distributions are in equilibrium with environmental conditions and that these species environment relationships can be used to estimate species responses to environmental changes ( Elith and Leathwick 2009 ). Relaxing this equilibrium assumption has been identified lately as key challenge in the field of species distribution modeling ( Dormann 2007 , Elith and Leathwick 2009 ). One way to begin addressing this challenge is to separate potential habitat estimation from the realized species distribution ( Soberon 2007 , Soberon and Nakamura 2009 ). Under non‐equilibrium conditions, habitat availability (areas in which environmental conditions allow species presence) and occupancy or abundance patterns differ ( Brotons et al. 2004 ). Recent studies have adopted hybrid SDM approaches combining correlative, phenomenological models with more mechanistic models, thereby allowing for the integration of critical ecological processes into SDMs ( Gallien et al. 2010 ). These new approaches purport to be more ecologically sound and achieve a better fit to the observed dynamics in species distributions ( Keith et al. 2008 , Elith and Leathwick 2009 , Willis et al. 2009 ). Such hybrid models explicitly incorporate ecological processes other than habitat selection (e.g. movement or competition), to better reflect actual species distribution patterns. Hybrid models allow us therefore to consider the dynamic processes of species distribution needed to model the effect of global change. As habitat connectivity and species dispersal capabilities are key factors determining a species use of newly appearing habitat and, hence, driving species responses to changes in habitat quality or availability it is important to include these processes in hybrid models ( Hanski and Ovaskainen 2000 ). However, dispersal has been explicitly incorporated in SDMs only in the last years ( Willis et al. 2009 ), and at mostly it has been used to trim model outputs using a rather static approach ( Peterson et al. 2001 , Allouche et al. 2008 ). These dynamical hybrid models will require however using datasets collected through time from validation purpose ( Keith et al. 2008 , Vallecillo et al. 2009 ). Here, we compare the ability of correlative and hybrid SDMs to quantify the distribution dynamics of an early succession bird, the ortolan bunting Emberiza hortulana , in a highly dynamic landscape driven by land abandonment and wildfire in Catalonia, north‐eastern Spain. We describe changes in this landscape using a landscape dynamics model that simulates fire impacts in the landscape and forest regeneration. Our hybrid model integrates the tracking of habitat suitability changes from a correlative model with a spatially‐explicit population model that incorporates dispersal, and assumes equilibrium in the initial state only. We investigate alternative scenarios of dispersal capability to generate predictions of distributions in the context of change in habitat suitability driven by fire disturbance. We finally validated correlative and hybrid model results using two independent data sets: one obtained from monitoring efforts deployed in recent burns; and the other obtained from a general‐purpose monitoring program, thus allowing models to be evaluated using a temporal perspective. The comparison of correlative and hybrid model performance allows us to test the degree to which species distribution dynamics (i.e. colonization patterns) are constrained by dispersal and/or habitat availability and thus evaluate the expected gain in model performance achieved by moving from purely correlative to hybrid SDMs in the case of species inhabiting highly‐dynamic landscapes. Methods Bird data The ortolan bunting is an open habitat, early‐succession bird species typically found in montane and subalpine meadows, open shrublands, as well as farmland, with a strong preference for recently burnt areas ( Brotons et al. 2008 , Menz et al. 2009 ). However, the analyses of large scale atlas data suggest that the use of recent burns may be constrained by dispersal thus making this species a good candidate to investigate species distributional changes in non‐equilibrium contexts ( Brotons et al. 2008 ). We compiled a range of large scale monitoring data relating to the distribution and spatial dynamics of the species in the region. First, comprehensive snapshot distribution data were obtained from the Catalan breeding bird atlas (CBBA) ( Estrada et al. 2004 ). The CBBA atlas sampled Catalonia (north‐eastern Spain) between 1999 and 2002 using a total of 3150 1 3 1 km plots in which presence/absence records were obtained from two one hour transects within each plot. Second, data on the species spatial dynamics were obtained from two bird monitoring programs from 2006 to 2009 in either unburnt (with relatively static habitat changes) or burnt (with more dynamic habitat change) sites. Information from unburnt sites was obtained from the Common Bird Monitoring Survey of Catalonia (CBS) ( Brotons et al. 2007 ). This general‐purpose monitoring program was initiated in 2002 and aims to track bird communities annually using 3 km transects divided in 500 m sections in which the number of observation of each species is noted at different distance bands. CBS transects are distributed within specific 10 3 10 km squares selected by volunteer observers throughout Catalonia. We compiled 3178 presence–absence records from CBS transects in 2006‐2009. Information from burnt sites was obtained from a similar monitoring method used to investigate bird species responses to recent wildfires in Catalonia (DINDIS) ( Zozaya et al. 2010 ). In this survey, started in 2006, all fires that occurred in Catalonia after year 2000 that are >50 ha are monitored yearly using a number of transects proportional to fire area. We compiled 1914 presence–absence records from DINDIS transects conducted between 2006 and 2009 on the 56 fires that occurred between 2000 and 2007. Environmental data and landscape dynamics Catalonia is a region with Mediterranean climate that is remarkably heterogeneous due to sharp climatic and geological gradients, including a range of landscapes from alpine habitats to coastal marshes, and from evergreen forests to steppes and agricultural mosaics. To describe habitat diversity relevant to our bird study model, we used geographical information system tools, to gather and combine information from available forestry, land cover and fire data ( Ministerio de Medio Ambiente 2006 ) to identify forested areas, their dominant species, shrubland and non‐forested open areas. The resulting forest map is composed of 50 3 50 m cells belonging to different forest (i.e. dominant species) and land use classes. Land abandonment and large wildfires appear to have driven regional landscape dynamics in recent years ( Díaz‐Delgado et al. 2004 , Gil‐Tena et al. 2009 ). We simulated temporal changes in landscape associated with these two processes using a spatially‐explicit landscape dynamic model implemented in the SELES (spatially explicit landscape event simulator) modeling platform ( Fall and Fall 2001 ). The landscape model includes a fire sub‐model and a shrub succession sub‐model. The fire sub‐model uses observed burnt areas between 2000 and 2009 directly and post‐fire vegetation changes were only modeled within these areas. Fire information was obtained from the historical sequence of fires available in Catalonia (Catalan Government, < http://mediambient.gencat.net/cat/el_medi/incendis/ >). Forest maturation was tracked using the years since the last fire. Post‐fire changes in forest composition (i.e. changes in dominant tree species or changes to shrub) and the probabilities for those transitions were obtained using bibliographic sources ( Rodrigo et al. 2004 ). We induced spatial autocorrelation in forest composition changes by coupling transitions of some cells with those of their neighborhood. Succession mimicking land abandonment impacts was only modeled for shrub (i.e. forests were assumed to be in a stable state). Shrub had a yearly probability of changing to forest after 20 yr. In case that a transition from shrub to forest occurred, the resulting forest type was chosen by using probabilities obtained from the distribution of mature forest types in the surroundings cells (150 m radius). The landscape simulation model was run 10 times for the period between years 2000 and 2009. Using the sequence of forest maps, produced by the landscape model, and additional topographic data, we then defined a set of landscape variables from which bird habitat suitability models could be derived (Supplementary material Appendix 1, Table A1). As the landscape model operates on a 50 3 50 m resolution grid, we averaged cell values for these landscape variables within a 1 3 1 km resolution grid to match the bird atlas and monitoring data. Correlative species models We developed three correlative SDMs of the distribution of ortolan buntings in the period 2000–2009 on a 1 3 1 km resolution grid ( Fig. 1 ). These models were based on different assumptions of species responses to environmental variables as follows. 1 Conceptual representations of the species distribution models used. The three correlative species distribution models; initial probability of occupancy (IPO); habitat suitability model (HS) and the combined habitat suitability model (CHS) were based statistical representations of habitat suitability, whereas the hybrid species distribution model included population dynamics and dispersal allowing habitat suitability to be tracked independently of occupancy. 1) Initial probability of occupancy (IPO) This model hypothesizes a static distribution throughout the studied period based on initial estimates of the species probability of occupancy based on the probability of occupancy for ortolan buntings in year 2000 ( Fig. 2 ), taken from ( Brotons et al. 2007 ). This correlative SDM included both environmental information as well as the spatial structure of bird distribution derived from the CBBA data (i.e. using autocovariates or other spatial predictors in the model) ( Augustin et al. 1996 , Lichstein et al. 2002 ). We assume that the IPO provides a good approximation for habitat suitability at year 2000. 2 Initial probability of occurrence (IPO) for ortolan bunting in Catalonia, taken from Brotons et al. (2007). 2) Pure environmental habitat suitability (HS) Given the hypothesis that habitat selection in the ortolan bunting is at least partially driven by changes in open habitats related to fire ( Menz et al. 2009 ), we used land use and forest variables describing burn areas and open‐habitat availability as explanatory factors (Supplementary material Appendix 1, Table A1) in a generalized linear model with binomial response. Model selection was based on AIC. To avoid strong departures from equilibrium in the training data, we calibrated the HS model using only CBBA data sampled within 10 km of a recorded presence. The highest ranked model included variables representing the percentage of young forest (<20 yr), and the percentage of unburnt and burnt shrubland in the 1 3 1 km cell (Supplementary material Appendix 1, Table A2). As shrubland and forest variables are updated by the landscape dynamics model, the HS derived an expected probability of species presence following changes in the landscape resulting from fire disturbance and succession processes. 3) Combined habitat suitability (CHS) In this third model, the IPO and HS models were combined to represent the effect of landscape changes while preserving information regarding the initial probability of occupancy. That is, we used IPO as an approximation to initial habitat suitability (CHS 0 =IPO) and used changes in HS to indicate increases or decreases of local habitat suitability. This leads to the following expression for such combined habitat suitability (CHS) for each cell at time t >0: CHS t = CHS t−1 + (HS t − HS t−1 ) where HS t is the value of the pure environmental habitat suitability model at time t . To arithmetically operate with the values of the IPO and HS, both models were expressed as probabilities. In the event that IPO=HS 0 , the environmental model would account for the species spatial structure and the equation above would simplify to CHS t =HS t . Hybrid species distribution modeling Hybridization of phenomenological with mechanistic models ( Gallien et al. 2010 ) facilitates the estimation of effective occupancy, defined as areas in which the species is likely to be present given both habitat suitability and ecological constraints on distribution limiting occupancy (e.g. dispersal limitation, competition) or favoring use of non‐favorable habitat (e.g. spillover effects). Here, we hybridized the combined habitat suitability (CHS) model with a spatially explicit population model. Assuming initial equilibrium, the hybrid model modifies an initial occupancy state at a given site (IPO) by tracking local changes in HS and the occurrence of the species in the neighborhood ( Akçakaya et al. 2005 ). Decreases in habitat quality will result in a progressive decrease in the individuals’ carrying capacity at that location and eventually resulting in local extinction. Increases in habitat quality in unoccupied sites (i.e. disturbance effects) lead to increases in the carrying capacity at that location, but these areas can only be colonized from neighboring occupied cells (i.e. population sources). Dispersal is, therefore, implemented as a key ecological constraint and leads to a possible delay in the occupancy of suitable habitat. When dispersal constraints are strong and habitat changes are fast, non‐equilibrium situations can arise rapidly. Although other mechanistic processes governing distribution, such as species interactions ( Willis et al. 2009 ), could be included, previous studies identify dispersal as the main constraint on ortolan bunting distribution ( Brotons et al. 2008 ). Like the spatially explicit fire‐landscape model, the spatially explicit population model was implemented in the SELES modeling platform ( Fall and Fall 2001 ). The population model is an age‐structured stochastic population model that tracks occupancy, abundance and carrying capacity for ortolan bunting males within each cell. The population model requires several inputs data and parameters which are described below. Initial occupancy and habitat suitability Initial occupancy and habitat suitability patterns are stochastically determined using the probabilities indicated in the correlative models. For any given cell i , a random uniform number between 0 and 1 is drawn and this number is compared to the IPO i value to determine initial cell occupancy. The same random number is compared to CHS i t to determine cell suitability at any given time t . Using the same random number ensures that: 1) there are not occupied cells in unsuitable habitats at time 0 (because CHS 0 =IPO); 2) cells are only allowed to change their suitability when changes in CHS occur (i.e. if CHS value does not change from t to t+ 1, the suitability state at t is the same as the state at t+ 1). Initial abundance and changes in carrying capacity Within the initially occupied cells the number of males in the cell was determined from a Poisson distribution with a mean based on a step function relating the initial probability of occurrence with density data from the breeding bird survey of Catalonia (CBS) from the 20022004 period ( Table 1 ; Brotons et al. 2007 ). The initial number male juveniles is then determined using a constant adult:juvenile ratio. The carrying capacity ( K ) at any time t was set using the same step function ( Table 1 ), but with CHS t as input instead of IPO. Whenever a previously suitable cell becomes unsuitable then K is set to zero, and if the cell was occupied then its population was forced to disperse. 1 Parameterization of the spatially explicit population model (SEPM) for ortolan bunting. Parameter Value (s) A) Initialization (from Brotons et al. 2007 ) (0.01–0.02): 6 IPO: abundance relationships (0.02–0.07): 8 (occurrence probability): (number of individuals) (0.08–0.22): 15 (same for CHS‐K relationships) (0.23–1): 18 B) Demographic parameters (from Steifetten and Dale 2006 ) Initial adult/juvenile proportion 0.60 Breeding probability N(0.65, 0.0045) truncated to avoid values outside [0,1] Offspring mortality probability 0.146 Juvenile mortality probability 0.67 Adult mortality probability 0.375 Probability distribution of a clutch of size j ( j =0, 1, …, 4) [0, 0.02, 0.14, 0.43, 0.41] Male proportion at birth 0.50 C) Dispersal parameters (following Dale et al. 2005 , 2006 ) Proportion of dispersal for adults 0.50 Proportion of dispersal for juveniles 1 Parameter α for negative exponential dispersal kernel D1: 1; D2: 0.5; D3: 0.355; D4: 0.1 Local population dynamics We designed and parameterized mortality and reproduction for ortolan bunting according to the literature ( Steifetten and Dale 2006 ) and local data. The first local breeding occurs at the end of the first year (in autumn), followed by local mortality (applied in winter). Reproduction and mortality events occur once per year. The probability an individual reproduced varied stochastically following a truncated normal distribution and, if reproduction occurred, the clutch size was determined from a multinomial distribution ( Table 1 ). For each newborn, the number of male offspring increased by one if the sibling is male and does not die after birth. The number of male juveniles for a given year is set to the male offspring, and the surviving juveniles from the last breeding event become adults. Different mortality rates were applied to adults and juveniles, which determined the number of individuals surviving the winter in each occupied cell. Note that this number may be bigger than the local carrying capacity. Dispersal Individuals disperse into neighboring cells with adequate habitat according to a dispersal kernel function ( Chapman et al. 2007 ). In most scenarios (see below), we modeled the distribution of dispersal distances (i.e. the dispersal kernel) as a negative exponential function of geographical distance ( d ): e ad , where the inverse of parameter α can be interpreted as the mean dispersal distance. To speed up computations, we truncated all kernel values below 0.001 to 0. The first dispersal event occurs one year after the simulation starts, and subsequent dispersal events are scheduled for every summer thereafter. According to the available information on the species dispersal capabilities ( Dale et al. 2005 , 2006 ), the proportion of individuals dispersing from each cell is 50% for adults and 100% for juveniles ( Table 1 ). However, if remaining adults still exceed the current carrying capacity, the additional individuals are forced to migrate. The location to where individuals can spread is limited to suitable cells that can be reached according to the dispersal kernel function and where the population is below the carrying capacity. We assume that habitat type does not affect dispersal ability. When there are no suitable locations available the disperser is assumed to die due to the unavailability of resources in the surroundings of the source cell. We included an additional parameter in our model to account for behavioral constraints in dispersal induced by conspecific attraction leading to biased dispersal to areas where the species is abundant ( Boulinier and Danchin 1997 , Keitt et al. 2002 ). We modeled conspecific attraction by multiplying the dispersal kernel of occupied cells by a constant factor (see below). Dispersal scenarios and model runs We ran the population model under eleven different dispersal scenarios: no dispersal (ND) – all dispersed individuals stay in the source cell, and those exceeding the carrying capacity of the cell are assumed to die; restricted dispersal (D1D5) – in these scenarios we tested three different values of the negative exponential parameter. In scenario D1 (α=1, 1/α= 1 km) dispersers can move up to 7 km (that is, with α=1 kernel values are lower than 0.001 for distances >7 km and we truncated those kernel values to 0), in D2 (α=0.5, 1/α=2 km) up to 13 km and in D3 (α=0.355, 1/α=2.81 km) up to 18 km. The distribution of dispersal distances in D3 reflects that observed in the ortolan bunting population in Norway ( Dale et al. 2005 ); long dispersal (D4D5) – in these two scenarios dispersers could move up to 50 km. In scenario in D4 we still used the negative exponential kernel (α=0.1, 1/α=10 km), while in D5 the distribution of dispersal distances was uniform between 0 and 50 km; treatments with conspecific attraction (D1CD5C) – we replicated the previous scenarios (D1D5) by adding an effect of conspecific attraction to dispersal. Specifically, we multiplied by five the value of the dispersal kernel of cells that were already occupied, whereas on suitable but unoccupied cells the dispersal kernel value was not modified. The five‐fold increase was chosen based on a sensitivity analysis that showed weak effects for small values. The spatially explicit population model was run for the period of 2000 to 2009 with 100 replicates for each combination of landscape simulation and dispersal scenario, and replicates were averaged to obtain the overall probability of occupancy. Model runs and validation Correlative and hybrid SDM predictions were validated using the two monitoring datasets. Observed occupancy in non‐burnt sites (3178 presence–absence records) obtained from the CBS monitoring program primarily corresponds to areas with small decreases or no change in habitat suitability. Observed occupancy data (1914 presence–absence records) obtained in burnt sites from the DINDIS project corresponds to areas affected by large wild fires ( Zozaya et al. 2010 ), where habitat quality is predicted to strongly increase and eventually lead to species colonization ( Menz et al. 2009 ). We used two different validation approaches to assess its predictive accuracy. In the first approach we calculated the agreement between the modeled probability of occurrence and the observed occurrence using area under the receiver operating characteristic curve (AUC) statistics ( Fielding and Bell 1997 ). In the second approach (Supplementary material Appendix 1, Table A3) we estimated the log‐likelihood of the model probabilities (in a binomial distribution) given observed occurrence data. We also estimated the coefficient of determination for the model (R 2 ) by comparing its log‐likelihood value with the log‐likelihood of a null model where the species prevalence is the only information available ( Nagelkerke 1991 ). Validation statistics were computed for two periods (20062007 and 20082009) taking within each period the maximum observed and modeled values among the two years. As time progresses and fire impacts accumulate in the landscape, we expect a progressive divergence in model predictive ability between correlative approaches and hybrid models including species dispersal constraints. Results Landscape, habitat suitability and occupancy dynamics Given that in our landscape habitat quality was generally decreasing (i.e. widespread land abandonment leading to a loss of open habitats), the total amount of suitable cells in the HS model, and hence in the CHS model, generally decreased ( Fig. 3 ). Only small increases were observed in years where the general decrease in availability was compensated for by the creation of new habitat due to fires. In the scenario in which a species could not disperse from the initially occupied areas (ND in Fig. 3 ), the total amount of occupied cells decreased following decreases in habitat suitability due to encroachment, forest succession and ageing. In the high dispersal scenario (D5 in Fig. 3 ), new habitats created by fire could be reached and, therefore, allowed the occupancy of new patches leading to progressive shift in distribution from areas initially occupied, but losing quality, to these post‐fire newly created habitats. In the scenarios in which dispersal occurred and individuals had the capability to move and find new habitat, occupancy patterns were between the patterns showed by the two extreme behaviors described above (D1 and D5 in Fig. 3 ). The higher the dispersal capability, the closer the match between real distribution patterns and habitat suitability changes will be. When dispersal is constrained, species distribution showed non‐equilibrium situations and increasingly differed from habitat suitability estimates (D1 and CHS in Fig. 3 ). Furthermore, our results indicated that conspecific attraction acted by effectively constraining dispersal and limiting occupancy to newly appearing habitat (e.g. compare D1 and D1C in Fig. 3 ). 3 Predicted habitat suitability and occupancy dynamics for the hybrid model under different dispersal scenarios. The thick solid line represents the number of sites with potential habitat in the region (i.e. sum of CHS values) in each year. The thin solid line represents the predicted number of occupied cells in each year under a non‐dispersal scenario (ND), which indicates losses of occupancy associated with decreases in habitat quality. The remaining dashed lines represent occupancy dynamics under different scenarios of dispersal capability, with and without conspecific attraction (scenarios D1C, D3C and D5C and scenarios D1, D3 and D5 respectively). Dispersal capability allows the species to partially compensate for the loss of occupied sites by making use of new suitable habitat associated with fire impact. Model validation results Among all comparisons, the HS model performed worse than the CHS correlative model ( Fig. 4a and Supplementary material Appendix 1, A3), indicating that considering habitat suitability changes starting from an accurate initial species distribution was more effective than a pure environmental niche model. In cases in which species data were gathered from the monitoring of non‐burnt sites, we detected weak differences among dispersal scenarios of the hybrid model and a good overall performance of the IPO model ( Fig. 4a , Supplementary material Appendix 1, A3). Moreover, although all the models were satisfactory in terms of agreement with observed data, the models including distribution dynamics performed worse than the scenario assuming a static distribution (IPO). Thus, changes affecting habitats already occupied in the initial conditions were very limited and hence initial occupancy predictions fit monitoring observations adequately. 4 AUC evaluation results of the hybrid model under different dispersal scenarios, using either the (a) common bird survey data (CBS) or (b) monitoring data from burnt habitats (DINDIS). Correlative models (IPO, HS and CHS) are indicated using black symbols. Dispersal scenarios are grouped into five categories: no dispersal (ND), restricted dispersal (D1D3), restricted dispersal with conspecific attraction (D1CD3C), long dispersal (D4D5) and long dispersal with conspecific attraction ( D4C D5C ). Strong differences between dispersal scenarios occurred for DINDIS evaluation data (20062007: F 4,105 =176.6, p<1 e‐15; 20082009: F 4,105 =1141.4, p<1 e‐15), but not for the CBS evaluation data (20062007: F 4,105 =3.4, p=0.0195; 20082009: F 4,105 =2.0, p=0.1208). However, when we used monitoring data from burnt habitats, the scenario predictions showed more marked differences both in the likelihood of the observed data and the predictive power to capture the ortolan bunting presence and abundance patterns after fire impact ( Fig. 4b , Supplementary material Appendix 1, A3). The scenarios that were most successful at predicting which burned areas were colonized by buntings included restricted dispersal capabilities (scenarios D1D3 and D1CD3C); both the no dispersal (ND) and long dispersal (D4/D5 and D4C/D5C) scenarios were less successful at prediction ( Fig. 5 ). Furthermore, in accordance with our prediction, we observed that differences between models increased in time, with outputs from the hybrid SDM with restricted dispersal capabilities offering the best fit to observed data ( Fig. 4b ). Thus, ortolan buntings appeared to have the ability to colonize newly appearing suitable habitats, but these capabilities were limited by dispersal constraints, leading to non‐equilibrium situations between species distribution and habitat suitability. 5 Probability of occurrence on fire patches in 2009 according to the spatially explicit population model for dispersal scenario D1 and landscape (run no. 3). Bird symbols indicate those fire patches where ortolan bunting was observed in monitoring data on burnt sites (DINDIS project). Discussion We have shown for the ortolan bunting that dynamic species distribution modeling based on the explicit integration of habitat suitability and effective occupancy estimation may successfully capture species distribution dynamics even when species distributions are not in equilibrium with habitat conditions. This is achieved by including key ecological processes that appear to limit effective occupancy of suitable habitat (i.e. dispersal constraints). Our hybrid modeling approach facilitates the inclusion of such processes by using initial species distribution information and a careful estimation of potential colonizer sources, and therefore dispersal fluxes, into new appearing habitat. New habitat arises as a function of explicit habitat selection hypotheses coupled with a landscape dynamics model. Although hardly surprising, our models suggest that when changes in landscape suitability are abrupt, a hybrid approach, including dispersal limitations, is superior to purely correlative dynamic approaches that project habitat suitability under the new environmental conditions. The benefit of hybrid dynamic SDMs over purely correlative SDMs may be especially apparent in the cases in which validation data is used to assess dynamic distribution models ( Araújo et al. 2005 ). Regardless, models should be able to reproduce observed dynamics in areas with sharp changes in habitat suitability and not only overall distribution patterns dominated by areas with more stable occurrence patterns. In this context, we encourage the validation of dynamic distribution models using temporally replicated biological data ( Araújo et al. 2005 , Willis et al. 2009 ). Projections of distribution change that are not validated with temporally replicated data should be used with caution and interpreted as a general optimistic hypothesis of potential distribution change ( Gregory et al. 2009 ). Such projections have been previously applied to the investigations of bird responses to changes in fire regime using a metapopulation modeling approach linked to a landscape simulation model similar to ours ( Akçakaya et al. 2005 ). However promising, without independent temporal validation of the models it may become difficult to assess the relevance of the model output in terms of the validity of the assumptions specifically implemented into dynamic SDMs. To date, hybrid modeling approaches have been used primarily as tools to assess the impact of climate change on the viability of populations and species distribution shifts ( Wintle et al. 2005 , Brook et al. 2009 , Willis et al. 2009 ). Future applications of hybrid species distribution models should allow the study of interactions between land use changes, climate change and changes in the associated perturbation regimes. However, the capture of ultimate processes allowing the track of changes in habitat quality for a given species and their match with adequate landscape models is probably not an easy task. In our case, the capability of our habitat suitability model (HS) to assess changes in habitat quality was constrained by the capability of the landscape model to describe landscape changes. These limitations have probably led the hybrid model to perform worse on stable areas than the initial probability of occupancy (IPO), causing the IPO to rank better in the validation exercises. Therefore, careful integration of initial distribution conditions, key ecological processes and changes in habitat suitability across space are, therefore, critical challenges ahead. Our results suggest that these factors will all be important and vary in importance in line with the ecological context considered. For the ortolan bunting, suitable habitat arises as a result of the alternative processes of land abandonment, leading to decreases in habitat availability, and fire, leading to the creation of new habitat patches, making assessment of the overall impacts on species distribution complex ( Brotons et al. 2008 ). Early succession species may respond in very different ways to such complex shifting mosaics depending of their ability to find newly available habitats ( Brotons et al. 2005 ). Dispersal has been shown to be a key parameter in the distribution dynamics of a species in a context of changing habitat availability ( Hanski and Ovaskainen 2003 ). Our models indicate that the ortolan bunting has a limited capability to colonize new habitats and is, therefore, only partially able to compensate for habitat loss in undisturbed habitats by colonizing new burnt sites. We suggest that limited dispersal capability could be driven in part by conspecific attraction that may restrict effective dispersal distances ( Serrano and Tella 2003 ). Model agreement with observed data was highly sensitive to changes in dispersal distance (directly, or induced by conspecific attraction). These results suggests that previous models incorporating dispersal in a static way may be flawed as dispersal cannot be separated from the distribution context in which it takes place ( Peterson et al. 2001 , Allouche et al. 2008 ). Given the limited dispersal capability of the ortolan bunting, and the widespread loss of open habitat, the spatial pattern of the new habitat is predicted to be an important driver of species persistence. Such a sinking–emerging island model, where new habitat suddenly appears and gradually vanishes in time, can be useful for a number of other species in a context of rapid environmental change. In such cases, restricted distribution ranges coupled with strong dispersal constraints will likely lead to non‐equilibrium situations and high variability in the distribution outcome depending on the spatial pattern of habitat change. Positive spatial autocorrelation in fire distribution will likely favor metapopulation persistence in the long term, whereas new distant fires may make new colonization difficult thereby decreasing chances of metapopulation maintenance. The prediction of species distribution dynamics is contingent upon estimates of dispersal capabilities. Our results suggest that, if adequate data exist, hybrid species distribution models could be used to estimate model parameters accounting for the dispersal ability better fitting the observed species spatial dynamics ( Smolik et al. 2010 ). This reverse modeling approach would be of similar nature to the rationale used to parameterized metapopulation models ( Lele et al. 1998 , Étienne et al. 2004 ). Given the accumulation of temporal biological data and adequate information on species distribution dynamics (e.g. monitoring schemes for a number of taxa; Yoccoz et al. 2001 ), this kind of approach may allow the estimation of dispersal capabilities for a large number of species. The time and spatial scale of the training data must also be appropriate to the rate of landscape dynamics and allow the capture of enough distribution changes since otherwise the dispersal ability of the species may be underestimated. To some extent the success in ecology of metapopulation theory and species distribution modeling ( Guisan and Zimmermann 2000 , Hanski and Ovaskainen 2000 , Guisan and Thuiller 2007 ) is rooted in the use of widely available biological data to parameterize spatial or temporal species distribution dynamics. Hybrid modeling approaches linked to the use of large scale temporally replicated biological data sets have the potential to become a key tool in setting the scenarios of species responses to global change impacts (Thuiller et al. 2008, Elith and Leathwick 2009 ). Acknowledgements We greatly acknowledge the task of all the contributors to the projects run by the Inst. Català d’Ornitologia, whose careful work has been indispensable for the development of the present study. We also thank Svein Dale for information on the ortolan bunting dispersal data and the GRAF (Catalan fire fighters brigade) for providing information on fire perimeters. Both Catherine Graham (the editor) and Hawthorne Beyer provided valuable comments: we are in their debt. This work has received financial support from the projects CGL2005‐2000031/BOS, CGL2008‐05506‐C02‐01/BOS and Consolider Montes CSD2008‐00040 granted by the Spanish Ministry of Education and Science (MEC). The study is also a contribution to the FP7‐226852 project SCALES. Partial funding was also provided by the Catalan Government grant SGR2009‐531 and a Beatriu de Pinós contract to MC (2009 BP‐B 00342).

Journal

Ecography – Wiley

Published: May 1, 2012

There are no references for this article.