Potential spread of the pine processionary moth in France: preliminary results from a simulation model and future challenges

Abstract

Context

Some forest insect pests are currently extending their range as a consequence of climate warming. However, in most cases, the evidence is mainly based on correlations and the underlying mechanism is not clearly known.

Aims

One of the most severe pests of pine forests in Europe, the pine processionary moth, Thaumetopoea pityocampa, is currently expanding its distribution as a result of climate warming and does not occupy entirely its potential habitat. A model describing its spread was developed to simulate its potential range in France under various climate change scenarios.

Methods

The spread model was divided into several sub-models to describe the growth, survival and dispersal of the species. The model was validated on the observed change of species distribution, its sensitivity was tested, and spread scenarios were simulated for the future.

Results

The model shows that climate warming initiated the species range expansion in France since the early 1990s. The spread is now limited by dispersal capability, but human-mediated dispersal could accelerate the range expansion.

Conclusion

Species range expansion is an indicator of climate change. However, time lags can appear due to limited dispersal capabilities, and human-mediated dispersal could create satellite colonies and artificially accelerate the spread.

1 Introduction

Since forests play an important role in carbon sequestration, assessing the effect of global change on forests is a key issue to better understand and predict future climate change (Bonan 2008). In addition to global change, biotic disturbance agents such as insects and pathogens can affect a forest’s health. Many of them are currently responding to climate change (Ayres and Lombardero 2000; Desprez-Loustau et al. 2006; Dukes et al. 2009; Netherer and Schopf 2010; Robinet and Roques 2010; Jactel et al. 2012), and their outbreaks could convert a part of the forests from a carbon sink to a carbon source (e.g. Kurz et al. 2008). This negative feedback is not the only one; economic losses resulting from lower wood production and health problems resulting from the outbreak of urticating species are others.

The pine processionary moth, Thaumetopoea pityocampa, is an insect which is currently extending its geographical distribution northwards and in higher elevation in France, and more generally in Europe (Battisti et al. 2005; Rosenzweig et al. 2007). This range expansion is caused by an increase of winter temperatures, resulting in better climate conditions for the larval feeding activity in winter and higher probabilities of winter survival (Battisti et al. 2005; Buffo et al. 2007; Robinet et al. 2007). This forest pest causes a decrease of pine growth (Jacquet et al. 2013) and an economic loss (Gatto et al. 2009). This species can spread not only by active flight of the moths but also by accidental transportation, generating long-distance jumps and the appearance of pioneer colonies far from the main distribution (Robinet et al. 2012). Since larvae can release urticating setae (Rodriguez-Mahillo et al. 2012), this forest pest has become an important urban pest on ornamental trees. Assessing its potential spread in response to climate warming and human-mediated dispersal could help us identify areas at risk, target more precisely the surveillance and the strategies to slow the spread, and also inform people about the sanitary risk just before pest arrival.

The potential spread of the pine processionary moth relies on several factors such as climate suitability, active dispersal capabilities, human-mediated dispersal and habitat distribution. The role of these factors is presented in more detail hereafter.

• 1. The first factor is climate suitability. Before the recent climate warming and the spread of the species, its geographical distribution was constrained by the combination of two climate variables: the monthly mean of minimum temperature in January and annual sunshine duration (Huchon and Démolin 1970). More recently, the daily mean of minimum temperature during the cold period was identified to be an indicator of the larval feeding activity and population survival (Buffo et al. 2007; Robinet et al. 2007). It is possible that, under changing conditions, a climate variable that could be neglected in the past (because of more important drivers) has now increasing importance.

• 2. The second factor involved in range expansion is the dispersal capability of the population. The observed spread rate was slightly above the known dispersal capabilities of the individuals. Female flight capability was thought to be very limited in the past (from a few hundred meters to 2 km; Démolin 1969), but higher flight capabilities were recently recorded in a flight mill experiment (up to 6 km for most of individuals, with a maximum of around 10 km for one individual; Robinet et al. 2012). Although this lab experiment may not reflect in total the flight capability in the field, these values are more consistent with the observed spread rate (5.6 km per year between 1992 and 2004; Battisti et al. 2005). Derived from the random walk theory, reaction–diffusion models have been used successfully to simulate the dispersal capabilities of many species (Shigesada and Kawasaki 1997). This type of models belongs to thin-tailed dispersal kernels which are particularly convenient for populations which exhibit a progressive loss of their genetic diversity along the colonization gradient and a genetic isolation by distance pattern (Klein et al. 2006). Since the genetic analysis of pine processionary moth populations shows such patterns in the main distribution (Rousselet et al. 2011), these models are particularly convenient to describe the dispersal capabilities of the moth.

• 3. The third factor playing a role in species range expansion is human-mediated dispersal. We therefore explore recent data on pioneer colonies, especially near Paris (Robinet et al. 2012), to determine the key factor and model these long-distance jumps.

• 4. In the potential spread of the pine processionary moth, another factor is involved: the distribution and density of host trees (mainly Pinus spp.). This host tree density is particularly important to calculate the carrying capacity and account for the effects of its heterogeneity over the territory.

Whilst some models have been developed independently to describe population dispersal (e.g. Roques et al. 2011) and population survival in relation with climate warming (e.g. Robinet et al. 2007), no model was currently available to describe simultaneously the potential growth and dispersal of the population at country scale under a changing climate for the pine processionary moth as for other species. Some models have been developed to describe short- and long-distance dispersal and population growth, for instance, to assess the potential spread of the horse chestnut leafminer, Cameraria ohridella (Gilbert et al. 2005); the potential spread of the pine wood nematode, Bursaphelenchus xylophilus, and pine wilt disease (e.g. Togashi and Shigesada 2006); and the potential spread of emerald ash borer, Agrilus planipennis (e.g., Muirhead et al. 2006). However, to explore the effects of climate change on species distribution, bioclimate models are generally used (e.g. Beaumont et al. 2007) instead of models describing the dispersal capabilities of populations. The objective of this study was to develop a novel type of spread model at country scale to assess more precisely the relative importance of climate warming and dispersal capabilities in the recent and future species range expansion.

2 Materials and methods

2.1 General framework

Hereafter, the unit of measure of the pine processionary moth density is the number of nests per pine to compare simulations to observations, but it is then converted into the number of nests per square kilometre in the model (given the number of pines per square kilometre).

We call the population density or distribution of year t the density or distribution of the generation hatching in year t − 1 and pupating in year t. To simplify the writing, the term “winter” is used hereafter to refer to the period from October of year t − 1 to March of year t (the period during which larvae generally encounter relatively cold temperatures that could affect their development and survival; Robinet et al. 2007). As the mechanism driving the extended diapause is not well known until now and the data collected so far do not allow modelling this mechanism, we do not consider the possibility of extended diapause over one or several years at the pupal stage.

The model described hereafter is applied over a grid of 1 × 1-km resolution covering France, and simulations were done with the statistical language R (R Development Core Team 2010). For each iteration, we applied in chronological order the following sub-models to the population density N _i,j (t) (number of nests per square kilometre), with (i,j) indicating the position on the grid and t the year of pupation. The list of parameters used in the model is given in Table 1; the list of input and output variables is given in Table 2.

Table 1 List of the parameters used in the model

Table 2 List of input and output variables

2.1.1 Larval survival as a function of minimum temperatures

Since a mismatch is observed between the historical model combining the minimum temperature in January and the annual sunshine duration (Huchon and Démolin 1970) and the current distribution (Fig. 1), we considered the recent indicator of the effects of climate warming: the mean of minimum temperature from October to March (Robinet et al. 2007). We used the survival function previously calibrated locally in the Paris Basin where latitudinal range expansion has been measured (Robinet 2006).

Fig. 1

Mismatch between potential distribution area defined more than 40 years ago and current distribution of the pine processionary moth in France. The black colour indicates the area where the species was supposed not to establish (called “exclusion area”, derived from Huchon and Démolin 1970), the grey colour indicates the area where the species is currently present (in winter 2010–2011), and the hatched area is the overlap between both areas. Grey dots with black point inside indicate the locations of pioneer colonies (all within the exclusion area)

$$ S(T)=\frac{1}{{1+\exp \left( {s\left( {{T_c}-T} \right)} \right)}} $$

(1)

with T the mean of minimum temperature from October to March (in degree Celsius), T _c the critical temperature (in degree Celsius) satisfying S(T _c) = 0.5, and s a sensitivity parameter. We used the same values for these parameters as for the Paris Basin: T _c = 3.3 °C and s = 1 (Robinet 2006). Since the nests are observed throughout winter, independently of their survival (some may contain larvae killed by harsh winter conditions), this survival function was applied to estimate the nest density after winter mortality (E _i,j (t)).

$$ {E_{i,j }}(t)={N_{i,j }}(t)S\left( {{T_{i,j }}(t)} \right) $$

(2)

This density provides a proxy of the population density at adult emergence in summer t.

2.1.2 Long-distance jumps

The analysis of pioneer colonies located far from the main species distribution revealed that they resulted from long-distance jumps due to human accidental transportation (Robinet et al. 2012). They are likely associated with the transportation and transplantation of large trees, such as ornamental trees generally planted in urban areas. Since this spread mechanism is a stochastic and relatively new event in time and space, modelling these long-distance jumps is challenging. To determine the factor that could explain the location of these pioneer colonies, we considered the host tree density (based on the IGN dataset—see details in “Growth model”) and the human population density (SEDAC 2012; Fig. 2). We calculated the 95 % confidence interval with a bootstrap method (with 1,000 samples) for the mean of these densities on the locations of the six pioneer colonies reported by Robinet et al. (2012) and on 100 points selected at random above the northern edge of the species distribution in 2010–2011. We also compared the survival rate in urban areas (defined as the area where the human population density exceeds 200 persons per square kilometre; Fig. 2a) and outside urban areas in France using the R function t.test (R Development Core Team 2010). To explore the effects of climate warming on the probabilities of establishment of transported individuals, we calculated the survival rate at these random locations for each year between 1980 and 2011 using the survival function S previously defined.

Fig. 2

Input variables. a Human population (persons per square kilometre). b Host tree density (trees per square kilometre)

Then, we developed a tentative model to simulate these long-distance jumps. The number of long-distance jumps was randomly chosen every year in the probability distribution Bin(n = 1, size = 1,000, p = 0.001) to simulate an average of one accidental transportation per year. This parametrisation is arbitrary since the parameters depend on various factors such as the number of infested trees transplanted each year coming from an infested area. Then, the destination of these long-distance jumps was randomly chosen in urban areas. We assumed that, for each long-distance jump, the equivalent of one nest was thus introduced in urban area independently of the distance to the main distribution; we corrected E _i,j (t) accordingly. This model, being only exploratory, was only used to draw scenarios for the future.

2.1.3 Local spread

Reaction–diffusion models have already been applied successfully to describe the spread of the pine processionary moth in the Paris Basin (Robinet 2006; Roques et al. 2011). We used the same model to simulate the potential spread in France. The initial condition was the population density at the time of emergence, E _i,j (t). The model is defined by the following formula:

$$ \frac{{\partial A}}{{\partial \tau }}\left( {x,y,\tau } \right)=D\cdot \Delta A\left( {x,y,\tau } \right)\;for\;\tau \in \left[ {0,1} \right] $$

(3)

$$ A\left( {{x_i},{y_j},\tau =0} \right)={E_{i,j }}(t) $$

(4)

$$ A\left( {{x_i},{y_j},\tau =1} \right)={A_{i,j }}(t) $$

(5)

with D the diffusion coefficient (in square kilometres per year), (x,y) the 2D continuous space, τ a subunit of time during adult dispersal, and A _i,j (t) the population density after active dispersal. This model was then discretized on the 1-km grid (Robinet 2006). Following a model calibrated to the observed spread south of the Paris Basin between 1972 and 2004, D = 2.9 km²/year (Robinet 2006). However, from Roques et al. (2011), D = 9.3 km²/year. This second model was calibrated on a shorter and more recent period (from 2007 to 2009) and did not consider climatic constraints. Both estimates were tested in the validation part, but only the first one was kept in the other parts.

2.1.4 Growth model

To account for a limitation of the population abundance (carrying capacity), no overlap of the generations and a delay effect in the plant–insect interaction, we considered the discrete form of the delayed Ricker model (e.g. Murray 2002) to describe the population temporal dynamics. The population density at time t + 1 was calculated as follows:

$$ {N_{i,j }}\left( {t+1} \right)={A_{i,j }}(t)\exp \left( {r\left( {1-\frac{{{E_{i,j }}\left( {t-1} \right)}}{{{K_{i,j }}}}} \right)} \right) $$

(6)

where r is the growth rate and K _i,j is the carrying capacity at position (i,j) on the grid. This carrying capacity depends on the pine density P _i,j , as follows:

$$ {K_{i,j }}=k{P_{i,j }}\Delta x\Delta y $$

(7)

with k the carrying capacity expressed as the number of nests per pine, P _i,j the pine density expressed as the number of pines per square kilometre, and ∆x∆y the area of the grid cells (which is equal to 1 km²). Based on field observations and preliminary data exploration (see Robinet 2006), r = 1 and k = 1 nest per pine were used in these simulations.

The GIS map of the Pinus spp. density was provided by IGN (French National Institute for Geographical and Forestry Information; http://www.ign.fr). Both Mediterranean pines (Pinus pinaster, Pinus pinea, Pinus halepensis) and mountainous pines (Pinus sylvestris, Pinus nigra, Pinus uncinata, Pinus cembra, Pinus mugo) were considered as potential host trees (Fig. 2b) despite differences in their susceptibility (Montoya 1981). Their density was measured in French forest stands covering more than 2.25 ha. The minimum host tree density recorded across forest stands in France was 11 pines/km². Isolated trees, especially ornamental trees in urban areas, are not surveyed at the national scale. However, they are widespread and could play an important role in the spread of the species in urban, agricultural or prairial areas.

Therefore, we assumed that a small density of pines is present outside the surveyed forest stands. From an exploratory analysis (data not shown), nearly all pines have a neighbour within a 3-km radius in an open-field region called “La Beauce”, where the area covered by pines (between Orléans and Paris) is negligible according to the IGN dataset (http://www.ifn.fr/spip/?lang=en). It means that the average density is at least one pine for 28 km². Therefore, we assumed that pine density was 0.04 tree/km² (Fig. 2b).

2.2 Simulations for model validation (1980–2011)

2.2.1 Initial conditions

The initial population distribution was the range observed in 1969–1979 (CTGREF-INRA 1980; Fig. 3). The population density in this range was supposed to reach the carrying capacity k. Since the population densities of two consecutive populations are needed to apply the model (see “Growth model”), we assumed that the initial population and the one before had the same density before winter. We then applied different mortality rates according to the temperature recorded in respective years.

Fig. 3

Distribution of the pine processionary moth in France downscaled to 1-km grid resolution for the model implementation: in 1969–1979 (in black; derived from CTGREF-INRA 1980), in winter 2005–2006 (in dark grey) and in winter 2010–2011 (in light grey). The two later distributions were derived from INRA surveys on a 8 × 8-km grid, with a maximum elevation at 1,400 m. Grey dots represent pioneer colonies discovered between 2003 and 2012 (derived from Robinet et al. 2012 and INRA latest surveys)

2.2.2 Climate dataset

We used the baseline temperature grid for the year 2000 (CCAFS 2012). Then, for each winter from 1978–1979 to 2010–2011, we calculated the difference between the mean of minimum temperature from October to March given by the baseline dataset and the mean temperature effectively recorded over ten weather stations homogeneously distributed in France (in the species range: Bourges, Merignac, Perpignan, Lyon, Rennes, Marignane, Orléans; outside the species range: Strasbourg, Metz, Beauvais; ECA&D 2012). We therefore used a temperature dataset mimicking the actual historical fluctuation in the model for the validation part (Fig. 4). In this way, the model simulations accounted for the very cold winter in 1985–1986, the warming trend (0.3 °C on average per decade between 1979 and 2011; +1 °C between 1970–1989 and 1990–2009) and the apparent temperature stagnation since the 2000s.

Fig. 4

Mean of minimum temperature from October to March in France between 1969 and 2011 (indicator of the larval feeding activity; Robinet et al. 2007). This temperature is based on the mean temperature in 2000 calculated over the grid covering France and the mean annual deviation recorded over ten weather stations. The grey line gives the yearly temperature and the black line gives the moving average over 10 years. Despite large fluctuations of the yearly temperature, the 10-year average temperature increased noticeably in the early 1990s and remained relatively stable since the early 2000s

2.2.3 Validation

To validate the model, we calculated the area under the curve (AUC) of the receiver operating characteristics plots (Fielding and Bell 1997). This measure is based on the sensitivity (the true positive rate) and specificity (the true negative rate) of the model independently of the density threshold chosen for predicting species presence. Then, we determine the maximum correct classification rate (CCR), defined as follows:

$$ CCR=\left( {TP+TN} \right)/\left( {TP+TN+FP+FN} \right) $$

(8)

with TP the number of true positives, TN the number of true negatives, FP the number of false positive and FN the number of false negatives. Presence/absence was derived from the simulated population density by assigning presence when the density was above the threshold density maximizing CCR. Since the northern edge of species distribution was rigorously surveyed at France scale during the winters 2005–2006 and 2010–2011 (Fig. 3; source: INRA Orleans), these distributions were chosen for validation. As previously mentioned, long-distance jumps were not simulated for the validation part; thus, pioneer colonies were not taken into account in this part of the study.

2.2.4 Sensitivity analysis

To assess the sensitivity of the model output to the parameters’ values, simulations were conducted from 1980 to 2011 with a change of each parameter ±25 % (one at a time). For each parameter value, we calculated the AUC, the CCR and the total invaded area for the year 2011. Due to large uncertainty about the value of the carrying capacity k, we also tested a higher value, k = 20 nests per pine.

2.2.5 Role of climate warming

To explore the role of recent climate warming in this northern range expansion, we compared the simulations of the spread model with the observed temperature fluctuations (like for the validation and sensitivity analysis parts) to simulations with historical temperatures corresponding to the average temperature between 1969 and 1979. This period corresponds to the historical distribution map at country scale (before the observation of the range expansion). To account for the parameters’ uncertainty, the values of the parameters were chosen at random in a normal distribution centred on the default values (see Table 1). The standard deviation of this normal distribution was chosen to provide around 90 % of the values within ±25 % of the default value (Table 1). The model was applied 100 times, and for each replicate simulation, we calculated the total invaded area and the y-coordinate of the most northern invaded cell in 2011. Presence/absence was derived from the population density given by the model simulation taking into account the density threshold found in the validation part. We used a Student’s test (t.test function in R; R Development Core Team 2010) to compare the mean of total invaded area and the mean of northern limit over the 100 replicate simulations between simulations with historical temperature (1969–1979) and simulations with observed temperature fluctuation (from 1980 to 2011).

2.3 Projections in the future (2011–2030)

2.3.1 Initial conditions

The initial population distribution for the projections in the future was the range observed in 2010–2011, including pioneer colonies far from the main distribution. The population density in this range was supposed to reach the carrying capacity k.

2.3.2 Climate change scenarios

Different climate change scenarios were taken into account to project this model in the future. We tested first an assumption of no more warming, taking into account the mean temperature between 1979 and 2011 (representing +0.9 °C compared to the baseline grid for the year 2000). Then, we tested different climate change scenarios for the 2020s. The potential spread was simulated under these scenarios until the year 2030. Although climate change scenarios are available for further time horizons, we focussed on near future only because the reliability of spread models generally deteriorates at each time step and their confidence level is relatively low over periods longer than a couple of decades (see Pitt et al. 2011). In terms of pest control and human health protection, the anticipation of pest arrival is needed mostly for the following years. In addition, beyond 10 or 20 years, no one can exclude possible adaptation of the species and different spread behaviour. We considered a widely used climate model (HADCM3 model; see for instance Beaumont et al. 2007) and tested emission scenarios supposed to cover a large range of socio-political–economic drivers of climate (A1B, UKMO-HADCM3 model; A2 and B2, HCCPR-HADCM3 models) (CCAFS 2012). Emission scenario B1 was not tested as it was not available when the datasets were downloaded. The mean of minimum temperature from October to March in France ranged from −8.9 to 10.4 °C, with an average of 3.0 °C in the previous scenario “no more warming”. Temperatures were on average 0.9 °C higher in scenario A1, 0.6 °C higher in scenario A2 and 0.8 °C higher in scenario B2. We simulated the potential spread in 2030 under these four scenarios with and without long-distance jumps. For these simulations, the values of the parameters were also chosen at random in a normal probability distribution centred on the default values with given standard deviations (see Table 1). The model was applied 100 times, and for each replicate simulation, we calculated the total invaded area invaded and the y-coordinate of the most northern invaded cell in 2030. Presence/absence was derived from the population density given by the model simulation taking into account the density threshold found in the validation part. We compared the mean invaded area and the mean of the most northern limit over the 100 replicate simulations for the eight scenarios (4 climate scenarios × 2 long-distance dispersal scenarios). We used a Student test (t.test function in R; R Development Core Team 2010) to compare these means and also calculated their 95 % confidence interval.

3 Results

3.1 Model validation (1980–2011)

The AUC = 0.84 and 0.75 and the CCR = 0.79 and 0.70, with D = 2.9 and 9.3 km²/year, respectively, for the year 2006; for the year 2011, AUC = 0.85 and 0.76 and CCR = 0.81 and 0.72. With the highest diffusion coefficient, the spread was slightly too rapid compared to the observed spread at the country scale.

3.2 Sensitivity analysis (1980–2011)

Changes in the value of the five parameters (Table 1) had various effects of the model output in 2011. The model was most sensitive to the threshold temperature (T _c), then to the growth rate (r) and, to a lower degree, to the carrying capacity (k) and the diffusion coefficient (D; Table 3). Although not very sensitive to the carrying capacity when testing the default value ±25 %, the model was relatively sensitive when testing a higher value. When considering a higher carrying capacity (k = 20 nests per pine), then AUC = 0.90, CCR = 0.86 and the total invaded area = 359,419 km².

Table 3 Sensitivity analysis (1980–2011): AUC, CCR and area invaded in 2011 when changing the parameters’ values

3.3 Role of climate warming (1980–2011)

Both the total invaded area and the northern limit of the species distribution were significantly higher in simulations with observed temperature fluctuation than in simulations with historical temperatures (1969–1979; total invaded area: t = 8.5619, df = 182.072, p < 0.001; northern limit: t = 3.7647, df = 197.704, p < 0.001; Fig. 5).

Fig. 5

Role of climate warming: mean of total invaded area (a) and mean of northern limit (y-coordinate in metric system “Lambert 2 étendu”) (b) in 2011 over the 100 replicate simulations with historical climate (1969–1979) and observed temperature fluctuation

3.4 Long-distance jumps

Human population density was significantly higher on the locations of the six pioneer colonies than on the 100 random points above the northern edge of the species distribution in 2010–2011 (CI_95% = 3,507–3,705 vs. 170–180). The mean host tree density over pioneer colonies was significantly below the density found on random points (CI_95% = 19–21 vs. 34–36). Over the six pioneer colonies, only one (17 %) was located in a forest stand with host trees (reported by IGN), which is slightly less than the probability to be there by chance. Therefore, human population density could be a key factor to define the potential locations for pioneer colonies, not host tree density.

Survival rate was slightly but significantly higher in urban areas in France (55 % in urban areas on average vs. 45 % outside; t = 79.98, df = 59134.2, p < 0.001), although the proportion of urban areas in France is relatively low (8 %). The potential survival rate on random points above the species range largely fluctuated between 1980 and 2011 (Fig. 6). However, three peaks of favourable conditions for establishing far from the main range clearly appeared in 1995, 2001 and 2007.

Fig. 6

Percentage of random points located in a favourable climate conditions (in case of long-distance jumps). The thin line represents a survival rate above 50 % and the thick line represents a survival rate above 75 %

3.5 Projections in the future (2011–2030)

As the total invaded area observed in 2011 was 339,121 km² and the northern limit was 2,415,025 m, the northern limit of the species range would continue to expand and move northwards even in case of temperature stabilization (Fig. 7). When taking into account climate change scenarios A1, A2 or B2, the results were very similar among them. The pine processionary moth could continue moving northwards slightly more rapidly than in case of no more warming. In addition, long-distance jumps increased considerably the species range whatever climate change scenario we considered.

Fig. 7

Projection in the future: mean of total invaded area (a) and mean of northern limit (y-coordinate in metric system “Lambert 2 étendu”) (b) in 2030 over the 100 replicate simulations with temperature stabilization (mean temperature over 1980–2011, Ø) and climate change scenarios A1, A2 and B2, with and without long-distance dispersal (LDD). Different letters indicate a significant difference (p < 0.05; Student’s t test) between the areas and the northern limits

4 Discussion

4.1 Potential spread

The model simulations clearly show the role of recent climate warming (between 1980 and 2011) in the species range expansion as the simulated spread was significantly higher when considering observed temperature fluctuation instead of historical temperatures (1969–1979; Fig. 5). The warming, particularly rapid in the early 1990s (Fig. 4), could probably explain the range expansion observed as soon as 1992–1996 south of the Paris Basin (Démolin et al. 1996; Goussard et al. 1999; Bouhot-Delduc 2005). In the future, climate change could still contribute to species spread. However, human-mediated dispersal could play a more important role in species range expansion and accelerate substantially the spread rate (Fig. 7). Therefore, the spread rate does not only depend on climatic constraints but also strongly depends on accidental transportation of the species. It also means that a large area beyond the species distribution edge is already favourable for species establishment (otherwise, transported populations could not establish and human-mediated dispersal would have no effect), and dispersal capability, either natural or human-assisted, is a crucial driver for the future. Even if accidental transportation of the pine processionary moth have probably occurred in the past, climate change increased the probability of these translocated populations to establish even when located far from the main distribution. Even if no pioneer colonies were discovered in the late 1990s (first peak of high survival; Fig. 6), they have been discovered in 2003, 2007 and 2008 (Robinet et al. 2012), in other words, just after the other peaks of high probability to survive above the species range (Fig. 6). Therefore, climate warming is involved not only in the continuous spread of the species but also in the possibility of stratified dispersal.

4.2 Climatic constraint

This study provides a preliminary assessment of the potential spread of the pine processionary moth in France. However, the model should be improved in the future. One of the most important parts to improve is the definition of the climatic constraint. Although the mean of minimum temperature from October to March was a good indicator of larval survival in the Paris Basin (Robinet 2006), this climate variable alone does not seem to explain correctly its presence in other parts of France, notably in mountainous areas (Massif-Central and the Alps). Determining more precisely the bioclimatic envelope of the pine processionary moth at large scale under a changing climate remains a challenge for the future. Although winter temperatures is the main driving factor of species range expansion (Battisti et al. 2005), other variables could constraint its bioclimatic envelope. For instance, the severe heat wave which occurred in France in 2003 is suspected to have killed locally a large part of the population (Bouhot-Delduc 2005; Piou et al. 2006). Indeed, the survival of winter populations could be affected by heat waves because temperatures higher than 36 °C have negative effects on young larvae (Santos et al. 2011). Understanding and modelling the impact of climatic anomalies in the species range expansion is very important in the future since the frequency of such extreme weather events is predicted to increase in the future (Beniston et al. 2007).

4.3 Need of a phenology model

Determining the climatic constraint at large scale is particularly difficult for the pine processionary moth because its phenology largely depends on the geographic region, as shown by Huchon and Démolin (1970). Therefore, the climatic constraint should be adapted to the phenology variability across France and, more generally, across Europe. Besides, the phenology scheme developed by these authors in the 1970s does not take into account changes in moth development induced by climate warming during the recent period. To assess the effects of climate change on species distribution, it is therefore important to understand how climate change will affect species phenology in the future. Consequently, the best approach would be to develop a phenology model accounting for climate change effects to determine the periods on which the climate variables should be considered for each life stage and to develop a fitness-based model. This type of model relies on both phenology and resistance to abiotic stresses and could contribute to better determine the species potential range (Chuine 2010).

4.4 Need of temperature projections at fine temporal resolution

Spread models are spatiotemporal models. Therefore, for each time step, the population density is calculated over a 2D space. Climate change scenarios generally provide an indication of the climate at several future time horizons. Temperature from climate change scenarios is generally not available for each year in the future. In addition to uncertainties associated with climate prediction for the future, these scenarios also introduce a bias in these spatiotemporal models as we assume a constant temperature between the current year and future time horizons. A weather generator, based on various climate change scenarios, could provide more realistic temperature fluctuation in the future. Considering fine temporal resolution datasets is all the more important as extreme weather conditions are predicted to occur more frequently and such conditions could eventually mitigate or enhance the effects of the mean increase of temperature.

4.5 Population growth and level

Another important mechanism to model is the growth of the population because the model is relatively sensitive to the growth rate. In addition, outbreak cycles are observed in the core area (area invaded for a long time), but not in the recently colonized area (Robinet 2006). Besides, the unit of measure of the population level is the number of nests per pine. This variable is a straightforward measure of the population on the field since nests are counted and then compared to the number of observed pines. However, this variable could introduce a bias to the real population level in terms of individuals (e.g. larvae). Indeed, a nest could contain various numbers of larvae (from a few tens to several hundreds). It is thus necessary to study the relationship between nest size, number of larvae inside and nest density.

4.6 Is the pine processionary a pure bio-indicator of climate warming?

Since the increase of winter temperature is the main factor which could explain the range expansion of the pine processionary moth, the northward shift of its northern edge was used as a bio-indicator of climate warming in France (ONERC 2012). However, climate warming velocity was higher than the species spread rate, and now the species does not occupy entirely the favourable area. A large part of France has become favourable for the species establishment, but the species has not yet spread over all this area (Robinet et al. 2007). Even in case of a warming stagnation (as we see since the 2000s; Fig. 4), the model suggests that the species could continue to spread. In addition, human-mediated dispersal could blur the signal and accelerate the spread. As a result, the range expansion of the pine processionary moth is an indicator of a more complex disturbance which includes climate change and human activities.

4.7 Advantages, limitations and perspectives of the model

The confidence level of simulations obviously depends on the validation of the hypotheses considered in the model and on the relevance of the functions used to describe the population dynamics. Therefore, simulations should not be considered as predictions but as scenarios that could happen under these assumptions. The originality of this model is to combine various sub-models that are generally developed independently (notably dispersal model, growth model and survival model in relation with climate warming). Although the model was sensitive to some parameters (see Table 3), it was useful to explore the role of climate change and human-mediated dispersal. This study revealed that integrating short- and long-distance dispersal, population growth and survival in a spatially explicit model applied to a heterogeneous landscape is possible and that it can provide a new tool to explore more precisely the role of different factors. This type of model can also be used to simulate the efficiency of control measures. Although not integrating the effects of climate change, Carrasco et al. (2010) combined natural dispersal, human-assisted long-distance dispersal and a phenology model to simulate the potential spread of the Western corn rootworm, Diabrotica virgifera virgifera. They explored the effects of applying buffer zones around newly infested areas and changing the rotation patterns of maize fields. Therefore, spread models like the one presented in this study can be used for testing different control strategies to contain the species or to slow its spread where the species is already present, or to assess the potential spread from hypothetical entry points where the species is not present.

4.8 Climate change effects at longer terms

As already mentioned, simulating the potential spread at the long term is problematic. In particular, we cannot exclude any potential adaptation of the moth populations, and especially of their biological cycle, to a new climate. Such adaptations to specific weather conditions prevailing in different biogeographic regions have been noticed for a long time (Huchon and Démolin 1970). However, the moth recently proved to be capable of large modifications in its biological cycle, such as was observed for the summer population in Portugal (Santos et al. 2011). Therefore, deeper studies are necessary to be developed in order to precise possible adaptations that could affect the moth life cycle. In addition, climate change will probably affect the natural distribution of host tree species at the long term (Cheaib et al. 2012). The way by which changes in host tree distribution could affect the geographical range of pine processionary moth is not clear because ornamental trees that can be used as hosts (pines, true cedars) are massively planted by humans not necessarily in the areas shown to be the most suitable for these trees. Another important unknown is the frequency of accidental transportation by humans in the future. Although this type of dispersal should be an important driver of the species distribution in the future, it is very difficult to know how this factor will change in the following decades.