MED-LMD Landscape dynamics modelling in the Meditteranean

Ortolan bunting hybrid SDM

Overview of the application

Ortolan bunting

The Ortolan Bunting, Emberiza hortulana (L.), is an open-habitat bird suffering large-scale declines across Europe due to general habitat loss and degradation as a consequence of farming intensification and homogenization of agricultural landscape. In contrast with overall European trends, recent surveys have demonstrated a recent remarkable population increase in Catalonia (NE of Spain) (Brotons et al. 2008). Although the Ortolan Bunting may be found in several habitat types, such as montane and subalpine meadows, open shrublands and farmlands, previous work has indicated the use of recently burnt areas in relatively high density, both in Catalonia and in other European countries. The analyses of large scale atlas data suggest that the use of recent burns may be constrained by dispersal thus making of this species a good candidate to investigate species distribution changes in contexts in non-equilibrium cases (Brotons et al. 2008).

Objective

In this application we compared the ability of traditional correlative and hybrid species distribution models to reproduce the distribution dynamics of the Ortolan Bunting in a highly dynamic landscape driven by land abandonment and wildfire in Catalonia.

Summary

Changes in landscape were obtained using the MEDFIRE landscape dynamics model. Correlative models can only project habitat suitability changes following landscape modification, hence assuming equilibrium with environmental conditions at all times. Instead, we used a hybrid species distribution model integrating the tracking of habitat suitability changes from a correlative model with a spatially-explicit population model including a dispersal module, and assuming equilibrium in the initial state only.

We ran our hybrid model under different scenarios of dispersal capability to generate different predictions of distribution change under a context of change in habitat suitability driven by fire disturbance. We then validated correlative and hybrid model results using two independent data sets obtained from independent monitoring efforts in recent burns and more stable habitats, thus allowing models to be evaluated using a temporal perspective. The comparison of correlative and hybrid model performance allowed us testing 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 moving from purely correlative to hybrid SDMs in the case of species inhabiting highly dynamic landscapes.

Landscape dynamic modeling details

We generated environmental dynamics related to fire and succession processes as follows. First, we took information from available forest maps in order to generate environmental layers to be used as predictors in habitat suitability models. Some of these layers were dynamically updated: the amount of forested areas, stable non-forested open areas, shrubland and non-forested open areas induced by fire. We used the MEDFIRE spatially explicit landscape simulation model (ver. 3, see Model archives) in order to update these dynamic environmental layers from 2000 to 2009. Since we knew the 2000-2009 historical fire sequence from available spatial databases, we used this information to determine areas impacted by fire disturbance. The succession sub-module was used to track post-fire transitions and forest maturation.

Species distribution modeling details

Overview of the hybrid SDM

We hybridized an habitat suitability model (HS) with a spatially-explicit population model (SEPM) parameterized for the Ortolan Bunting. Assuming initial equilibrium, the hybrid model modifies an initial probability of occupancy (IPO) by tracking local changes in HS and the distribution of the species in the neighborhood. Decreases in habitat quality will result in a progressive decrease in the individuals’ carrying capacity of the location and eventually lead to local extinctions. Increases in habitat quality in unoccupied sites (i.e. disturbance effects) lead to increases in the carrying capacity of the location, but those areas can only be colonized from neighboring occupied cells (i.e. population sources). Dispersal is therefore the key ecological constraint implemented and leads to variability in the potential of occupying suitable habitat. In cases in which dispersal constraints are strong and habitat changes fast, non-equilibrium situations can rapidly arise.

Initial Probability of Occupancy (IPO)

There are multiple ways in which environmental changes affect species distribution but one of the most apparent is by inducing changes in the quality and availability of resources required. Given the hypothesis that habitat selection in the Ortolan Bunting is driven by changes in open habitats related to fire, we used a set of land use and forest variables describing fire impact and open habitat availability, as well as broad climatic variables, as explanatory factors in a generalized linear model with binomial response. In order to avoid departures from equilibrium in the training data, we calibrated this PEHS model using only data sampled not further away than 10 km from a recorded presence. Starting from a model with all the landscape variables, we conducted an AIC-based backward selection (i.e. at each step we removed the variable leading to a larger decrease in AIC until all removals produced higher AIC values). The predictor of the final model included three variables: the percentage of young forest (less than 20 years), as well as the percentage of unburnt and burnt shrubland. Since shrubland and forest variables are updated within the landscape dynamics model, the PEHS derived an expected probability of species presence following changes in the landscape introduced mostly by fire disturbance and succession processes. We used a simple SELES model to translate changes in the landscape layers into changes of PEHS values (see HSI ver.2 in Application archives).

Habitat suitability (HS)

We combined IPO with changes in PEHS to track the effect of landscape dynamics on habitat suitability. That is, we used IPO as an approximation to initial habitat suitability (HS₀ = IPO) and use changes in PEHS to indicate increases or de creases of HS for a given species

environmental relationship hypothesis. This leads to the following expression for habitat suitability (HS) at time t > 0:

HS_t = HS_t_-1 + (PEHS_t – PEHS_t_-1)

where PEHS_t is the value of the pure-environmental habitat suitability model at time t. In order to arithmetically operate with the values of the IPO and PEHS, both models have to be based on probabilities. In the unlikely situation where IPO = PEHS₀, the pure environmental model would account for the whole species spatial structure and the equation above simplifies to HS_t = PEHS_t. We calculated HS from IPO and PEHS inside the SELES spatially-explicit population model described below.

Spatially-explicit population model

A stochastic population model was implemented in SELES modeling platform (see SEPM ver.5 in Application archives). The model tracks ccupancy and abundance at the cell level in terms of number of individuals and potential population value (i.e. carrying capacity).

Initial occupancy and habitat suitability

Initial occupancy and habitat suitability changes are stochastically determined. For any given cell i, a random uniform number is drawn. This number is compared to the IPO_i value in order to determine initial cell occupancy. The same random number is compared to HS_it in order 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; (2) cells are only allowed to change their suitability when changes in HS occur; and (3) the expected number of suitable cells at any time t is equal to the sum of HS_it values.

Initial abundance and changes in carrying capacity

Once a cell is set as occupied, the number of males in the cell is determined using a Poisson distribution. In the case of ortolan buntings, we set Poisson’s lambda using a step function relating the initial probability of occurrence (IPO) with available density data from the breeding bird survey of Catalonia (CBS) from the 2002-2004 period (Brotons et al., 2007). The initial number of adults and juveniles is then determined using a constant adult/juvenile proportion. The carrying capacity (K) at any time t is set using the same occurrence-abundance relationship, but with HSt as input. Whenever a previously suitable cell becomes unsuitable then K is set to zero, and if the cell was occupied then its individuals are forced to leave by dispersal. On the other hand, if a cell becomes newly suitable K is set as described.

Local population dynamics

We designed and parameterized mortality and reproduction for ortolan buntings according to the literature (Steifetten et al., 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). Successive breeding and mortality events occur once a year. The probability for an individual to reproduce varies stochastically following a normal distribution and, if reproduction occurs, a clutch size is picked from a multinomial distribution. For each newborn, the number of male offspring increases by one if the sibling is male and does not die after birth. The number of juveniles for a given year is set to the male offspring, and the surviving juveniles from the last breeding event become adults. The probability of death is applied differently for adults and juveniles and sets the number of individuals finally surviving winter in each occupied cell. Note that this number may be bigger than the local carrying capacity.

Dispersal

By modeling abundance within a cell, we were able to estimate a potential colonizer flux emitted by each occupied cell. Individuals disperse into neighboring cells with adequate habitat according to a dispersal kernel function, for which several options are possible. In most scenarios, we modeled the distribution of dispersal distances (i.e. the dispersal kernel) as a negative exponential function of geographical distance (d):

We considered a model of the probability of occupancy for ortolan buntings in year 2000 at 1x1 km resolution, taken from Brotons et al. (2007). This is an accurate correlative SDM that included both environmental information as well as the spatial structure of the distribution derived from real data (i.e using autocovariates or other spatial predictors in the model). If we assume initial equilibrium (i.e. species distribution determined uniquely by environmental variables), the IPO may also provide a good approximation for habitat suitability at time zero. Because it was built using the spatial structure of real (observed) data, the IPO model cannot be dynamically projected without including the factors that account for the spatial structure of the distribution. We used this model as a simple hypothesis of a static species distribution insensitive to environmental changesHabitat suitability modelling

Pure-environmental habitat suitability (PEHS)

Model evaluation

, where the inverse of parameter a can be interpreted as the mean dispersal distance. In order to speed up computations, we truncated all kernel values below 0.001 to 0. The first dispersal event occurs after a year, and subsequent dispersal events occur every summer thereafter. According to the available information on the species dispersal capabilities (Dale et al., 2005, Dale et al., 2006), the proportion of individuals dispersing from each cell is 50% for adults and 100% for juveniles. 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. When there are no suitable locations available the disperser is forced 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. We modeled conspecific attraction by multiplying the dispersal kernel of occupied cells by a constant factor.

Model execution and validation

Model scenarios

We ran the population model under eleven different dispersal scenarios:

- No dispersal (ND) – All dispersed individuals stay at the same cell they originate, and those exceeding the carrying capacity of the cell are forced to die.
- Restricted dispersal (D1-D5) – In these scenarios we tested three different values of the negative exponential parameter. Due to the truncation of low kernel values, in scenario D1 (a = 1, 1/a = 1 km) dispersers can move up to 7km, in D2 (a = 0.5, 1/a = 2 km) up to 13km and in D3 (a = 0.355, 1/a = 2.81 km) up to 18km. The distribution of dispersal distances in D3 tries to mimic the observed in the Ortolan population in Norway (Dale et al., 2005);
- Long dispersal (D4-D5) – In these two scenarios dispersers could move up to 50km. In scenario in D4 we still used the negative exponential kernel (a = 0.1, 1/a = 10 km), while in D5 the distribution of dispersal distances was uniform between 0 and 50 km;
- Treatments with conspecific attraction (D1C-D5C) – The previous five scenarios including some dispersal were replicated by adding an effect of conspecific attraction, which consisted in multiplying by five the dispersal kernel of already occupied cells with respect to that of suitable but unoccupied cells.

Model runs and evaluation

The spatially-explicit population model was run for the period of 2000 to 2009 having 100 replicates for each combination of landscape run and dispersal scenario, and runs were averaged to obtain the probability of occupancy. Model results were evaluated using the two temporal monitoring datasets. Observed occupancy in non-burnt sites obtained from the CBS monitoring program mostly corresponds to areas with little (i.e. decreasing) or no change in habitat suitability. Observed occupancy data obtained in burnt sites from the DINDIS project, in which bird colonization patterns were analyzed in large wild fires in Catalonia (Zozaya et al. in press), allowed a fire-colonization evaluation. In these burnt areas habitat quality is predicted to strongly increase and eventually lead to species colonization.

We evaluated each model’s results using the agreement between the modeled probability of occurrence and the observed occurrence using area under the receiver-operating curve (AUC) statistics. Evaluation statistics were computed for two periods (2006-2007 and 2008-2009, taking within each period the maximum observed and modeled values among the two years) and averaged over all landscape simulation runs.

Model results

Habitat suitability 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 PEHS model, and hence in the HS model, generally decreased; only small increases in years where the general decrease in availability was compensated by creation of new habitat due to fires. In a scenario in which a species cannot disperse from the initially occupied areas (e.g. ND), the total amount of occupied cells is expected to decrease following decreases in habitat suitability due to encroachment and forest succession and aging. In a hypothetical high dispersal scenario (e.g. D5), new habitats created by fire would be always reached and therefore allow 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 those cases in which dispersal does occur and individuals have the capability to move and find new habitat, occupancy patterns would fall between the patterns showed by the two extreme behaviors described above (e.g. D1 or D5). The higher the dispersal capability of the species, the closer the match between real distribution patterns and habitat suitability changes will be. In cases in which dispersal is constraining species distribution, the latter will show non-equilibrium situations and will increasingly differ from habitat suitability estimates. In addition, our results indicate that conspecific attraction will act as effectively constraining dispersal by limiting occupancy to newly appearing habitat (e.g. compare D1 and D1C).

In all comparisons the PEHS model performed worse than the HS correlative model, which indicates that considering habitat suitability changes starting from an accurate initial 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. Moreover, whereas all the models were satisfactory in terms of agreement with observed data, any of the models involving distribution dynamics performed worse than the scenario assuming a static distribution (IPO). Thus, changes affecting stable habitats already occupied in the initial conditions are weak and generally fit the data available. However, when we used monitoring data from burnt habitats, the scenario predictions showed marked differences both by the likelihood of the observed data and predictive power to capture the bunting presence and abundance patterns after fire impact. Both the no dispersal (ND) and long dispersal (D4/D5 and D4C/D5C) scenarios showed low power to predict fire areas that were effectively colonized by the species. Restricted dispersal capabilities (scenarios D1-D3 or D1C-D3C) were better at predicting species occurrences at new fires. In short, ortolan buntings do have the ability to colonize newly appearing suitable habitats, but these capabilities are limited by dispersal constraints and this leads to non-equilibrium situations between the distribution and habitat suitability of the species.

DINDIS

CBS

References

- Brotons L, Herrando S, Pla M (2007) Updating bird species distribution at large spatial scales: applications of habitat modelling to data from long-term monitoring programs. Diversity and Distributions, 13, 276-288.
- Brotons L, Herrando S, Pons P (2008) Wildfires and the expansion of threatened farmland birds: the ortolan bunting Emberiza hortulana in Mediterranean landscapes. Journal of Applied Ecology, 45, 1059-1066.
- Dale S, Lunde A, Steifetten O (2005) Longer breeding dispersal than natal dispersal in the ortolan bunting. Behavioral Ecology, 16, 20-24.
- Dale S, Steifetten O, Osiejuk TS, Losak K, Cygan JP (2006) How do birds search for breeding areas at the landscape level? Interpatch movements of male ortolan buntings. Ecography, 29, 886-898.
- Steifetten O, Dale S (2006) Viability of an endangered population of ortolan buntings: The effect of a skewed operational sex ratio. Biological Conservation, 132, 88-97.
- Zozaya, E., Brotons, L., Herrando, S., Pons, P., Rost, J. & Clavero, M. (in press) Monitoring spatial and temporal dynamics of bird communities in Mediterranean landscapes affected by large wildfires. Ardeola.

Page updated

Google Sites

Report abuse