Recorded and potential distributions on the Iberian peninsula of species of Lepidoptera listed in the habitats directive

Abstract. Using data on the known Iberian distributions of 10 species of Lepidoptera listed in the Habitats Directive referenced to the 10 × 10 km UTM grid, we determined their potential distributions and their relationships with selected bioclimatic factors associated with mean temperature and precipitation using Ecological Niche Factor Analysis (ENFA). Scores for Specialization and Marginality were determined in order to evaluate the relationships between the predictions of the model and climatic factors. The number of squares on the Iberian Peninsula in which the species are recorded and those squares predicted to be favourable for these species were determined if they matched the network of Protected Natural Areas. This suggested that a further eight 10 × 10 km squares should be included in Protected Natural Areas. The results also indicate that climate determines the distributions of most of the species. Although overall there is a close association between the observed and predicted distributions, the less thoroughly documented geographic ranges (i.e. those of the moth species) depart from this pattern.


INTRODUCTION
The Habitats Directive (HD, also known as Council Directive 92/43/EEC on the conservation of natural habitats and of wild fauna and flora) is one of the most important tools for the conservation of biodiversity in the European Union, which aims to preserve biodiversity by the conservation of natural habitats of fauna and flora.The designation of these habitats is strongly influenced by the biology and the distribution of the Species of Community Interest and Priority Species, both listed in Annex II of the HD.These species receive special attention and the data on their geographic distributions are being currently updated.The Annexes of the Habitats Directive include several orders of European animals and plants.Here we focus on Lepidoptera (moths and butterflies) within a specific territory, that of the Iberian Peninsula.Ten species of Iberian Lepidoptera are included in Annexes II, IV and V of the Habitats Directive and one of them is an Iberian endemic species: Polyommatus golgus (Hübner, 1813).
Lepidoptera are characterized by a huge diversity in terms of described species (more than 150,000 species within 128 families: van Nieukerken et al., 2011).Although the chorology of the butterflies (superfamily Papilionoidea) on the Iberian Peninsula is reasonably well studied (Romo et al., 2006), there is evidence that relevant areas remain insufficiently sampled (Romo & García-Barros, 2005) and for the faunistics of non-papilionoid Lepidoptera (moths in a wide sense) it is even worse.Given the present low investment in this kind of research, alternative techniques, such as distribution modelling, may be crucial tools for conservation purposes.Predictive models are now widely used in conservation biology (systematic conserva-tion planning), although decisions about conservation are currently mainly based on subjective criteria, such as personal experiences or expert judgment (Cook et al., 2010).Few applied conservation practices are based on results of using systematic conservation planning, despite the objectiveness of this method (Cowling et al., 2003).There are studies on the problems of modelling, such as the lack of certainty (Addison et al., 2013).A good approach would be to combine the objective results of models with expert knowledge, which is always to some extent subjective.The number of models available is overwhelming and the predictions vary depending on the model selected.Thus, it is essential that there is a good quality and refined database on the presence of the species, that environmental characteristics are predictable and the most suitable model is used (Lobo, 2008).Models do not have to be complex to be reliable.In fact, simple models with a few explanatory variables (between one and three) are more robust and generally preferred by experts (Van Couwenberghe et al., 2013).The available modelling methods include regression, machine learning or maximum entropy, among others (Elith et al., 2006;Elith & Leathwick, 2009;Tarkesh & Jetschke, 2012).We used Ecological Niche Factor Analysis, which is based on presence-only and compares the environmental values of the squares where the species has been recorded with those of the area studied (Hirzel et al., 2002).In this way, the most relevant variables are selected for the species and some uncorrelated factors that define Marginality (the distance between the optimal climatic conditions for the species and the average climatic conditions in the area studied) and Specialization (describing how specialized the species is based on the ratio between Lycaenidae); Eriogaster catax (Linnaeus, 1758) (family Lasiocampidae); Actias isabellae (Graells, 1849) (family Saturniidae); and Proserpinus proserpina (Pallas, 1772) (family Sphingidae).Three further species included in the annexes and recorded in the past from the study area were discarded: the Nymphalid Coe nonympha oedippus (Fabricius, 1787), Erebid Euplagia quadri punctaria (Poda 1761) and Lycaenid Lycaena helle (Denis & Schiffermüiller, 1775).C. oedippus does not occur on the Iberian Peninsula as far as we know as previous records are misidentifications (Fernández, 1991;Vives Moreno, 1994;García-Barros et al., 2013).E. quadripunctaria is found abundantly throughout the Iberian Peninsula (Ylla et al., 2010) and the subspecies that is protected is endemic to the island of Rhodes (Greece): Eupla gia quadripunctaria rhodesensis (Legakis, 1996).L. helle was recorded from just two 10 × 10 km UTM squares more than 25 years ago and therefore it seems necessary to confirm its occurrence in the study area (García-Barros et al., 2013).
The data on distributions came from several databases of presences (García-Barros et al., 2004) updated to 2010 and from a specific database created for the project on updating Habitats Directive species distributions (VV.AA., 2010).In total they provided 4,687 distribution records.

Climatic variables
We used 19 of the bioclimatic variables (BIOCLIM, Table 1) in Worldclim database (http://www.worldclim.org/bioclim) that represent gradients in mean temperature and precipitation, both seasonal or annual, or limiting climatic factors (Hijmans et al., 2005).These variables are the result of the collection of monthly measures of different variables related to temperature and precipitation, recorded at a large number of meteorological stations spread around the world, between the years 1950 and 2000.The resolution of five arc-minutes was adopted since it is the most similar to that of our species presence data (squares of 10 × 10 km).
Since our data consist of presences and several of the variables involved may be highly correlated with one another, the most relevant variables were chosen by using Ecological Niche Factor Analysis (ENFA) (Hirzel et al., 2002) implemented in Biomapper v. 4.0 (Hirzel et al., 2005).We consider the ecological niche (Hutchinson, 1957) as the subset of squares in the study area where the species have a reasonable probability of occurrence.This multivariate niche can be quantified in terms of both Marginality and Specialization index.From this point of view an ENFA is analogous to a Principal Component Analysis (PCA), in that it computes a group of uncorrelated factors with ecological meaning, and compares the climatic data at the localities of known presence of the species with the existing conditions in the study area.Marginality values obtained from the ENFA give an indication of the distance between the mean optimal value detected for the species for each factor selected and the average climatic conditions in the area studied, while Specialization values are the ratio between the climatic variability in the study area with that existing in the squares where the species selected is present.These values indicate the tolerance of each species to other secondary environmental gradients in the study area.Positive Marginality values indicate that the species studied prefers localities with values higher than the mean value for the whole area considered (Hirzel et al., 2002).The sign of the Specialization values has no biological meaning, although the variables with highest values are those most important in restricting a species geographical distribution.
The factors were retained after comparing their eigenvalues with those determined based on a broken-stick distribution (Hirzel et al., 2002).The variables selected were those with the high-climatic variability in the study area and that of in the area where the species occurs) determined.This method has a solid conceptual basis (Calenge & Basille, 2008) and has been successfully used in many other studies (Lobo et al., 2010;Chefaoui et al., 2011).
The existence of many species of Lepidoptera is under threat.The ranges of several species have contracted over the last few decades due to human activities, including climate change (Warren et al., 2001;Hill et al., 2002;Thomas et al., 2004;Stefanescu et al., 2011).In order to mitigate the effects of human pressure, a number of protected areas were created in Spain and Portugal (Protected Natural Areas, PNA).In these areas priority is given to the maintenance or restoration of a favourable conservation status for the habitats and species on the Iberian Peninsula.For some groups such as vertebrates (Rey Benayas & de la Montaña, 2003) the PNA does not include all the species or areas of high diversity value.On the other hand, most species of Iberian butterflies occur within the PNA (Romo et al., 2007).As the PNA are probably not equally adequate for all the animal groups, it is worth exploring their suitability for well known priority butterfly and moth species in the Iberian Peninsula, which has not been done yet.Of the species of moths listed in the DH only Actias isabel lae (Graells, 1849) has previously been analyzed from a biogeographic point of view in the study area (Chefaoui & Lobo, 2007).
For the reasons discussed above (the broad knowledge of their present distributions, unequal sampling effort, gradual reduction in their ranges), Lepidoptera are good candidates for species distribution modelling (Fleishman et al., 2001;Scheingross, 2007;Newbold et al., 2009).Therefore, it seems appropriate to model the distribution of the species of Lepidoptera listed in the Habitats Directive and identify potential areas with suitable habitats for the priority species included in the Directive.
The aims of this work, therefore, consisted in (1) modelling the distributions of species of Lepidoptera included in the Habitats Directive, (2) predicting possible new areas of species occurrence and (3) comparing these predictions with the network of Protected Natural Areas.This will enable us to suggest further areas for inclusion in PNA since they have favourable habitats for priority species of Lepidoptera living on the Iberian Peninsula.

Taxonomic group and study area
The study area is the Iberian Peninsula, using cells of 10 × 10 km based on the Universal Transverse Mercator (UTM) projection system for the representation and geocoding of the distributions of species.

Modelling
Maps of the potential distributions of the species were obtained using the occurrence data for each species of Lepidoptera as the dependent variable and climatic variables selected by ENFA as explanatory ones.The models were then fitted using the BIO-CLIM (Busby, 1991) program implemented in Diva-Gis v.7.1.7 (Hijmans et al., 2012).

Conservation status
Since all the species considered in this study are protected at an European level, using ArcMap 9.3 (ESRI, 2009) we seek matches with the Iberian PNA and (1) the number of 10 × 10 km squares in which the species are currently present and also (2) with those squares that the model indicates are potentially suitable for the species.In the first case a density map was produced by overlapping the species presence on the Iberian Peninsula to determine whether the known distribution of these priority species matches the priority areas for conservation.In the second case, mismatches are proposed as areas for conservation as the conditions in these areas are favourable for the species.The data on the protected areas in Spain came from Ministerio de Agricultura, Alimentación y Medioambiente (http://www.magrama.gob.es/es/) and Europarc (www.europarc-es.org)for 2013 and 2011, respectively.The data for Portugal were obtained from ICN (Instituto da Conservação da Natureza, http://www.icn.pt/) and Garcia-Pereira et al. (2003).
For practical reasons we assumed that a cell (UTM square) in which a species is present or predicted to host the species, was effectively protected (i.e., inside a protected area) if more than 15% of it was in the protected area.This resulted in a total of 1,210 UTM squares being included in PNA.

Actual distribution and density of species
Four of the species studied are comparatively widespread: A. isabellae, P. proserpina, E. aurinia and P. apol lo, the last two with the highest number of geographical records: 1,121 and 2,133 respectively.The remaining six species have restricted distributions: E. catax, L. achine, P. arion, P. golgus, P. mnemosyne and P. nausithous (Figs 1, 2).
The density map show that the squares with highest numbers of species are located in the Pyrenees, where there are two squares with six of the ten species listed in the Habitats Directive (Fig. 3A).These two squares (30TXN83 and 30TXN84) are located in the Valle de Hecho, Huesca province (northern Spain) and each hosts E. aurinia, E. catax, A. isabellae, P. apollo, P. arion and P. mnemosyne.

Bioclimatic factors
The ENFA analysis selected 12 of the 19 bioclimatic variables associated with precipitation and temperature, which accounted for between 85 and 99% of the distribu-tions of target species (Table 2); in all cases the ENFA analysis selected three uncorrelated variables for each species that were capable for explaining their distributions.
The variable selected in most cases was "Maximum Temperature of Warmest Month" (MATWM), which is significantly associated with the distribution of six of the ten species.This variable was followed by "Annual Precipitation" (APREC), the "Mean Temperature of Coldest Quarter" (MTCQ), the "Precipitation of Driest Month" (PDM), the "Temperature Annual Range" (TAR) and "Temperature Seasonality" (TS).All these variables appear independently three times as important factors associated with the distributions of the species.Variables APREC, MTCQ and TAR together explain the distributions of two of the species listed in the Habitats Directive: P. golgus and P. nausithous.
Marginality values indicate that the optimum climates for most of the species differ from the average climatic conditions of the peninsular area (Marginality > 1.5) and the climatic optimums of only three species (E.catax, P. golgus and P. nausithous) are close to the average climatic conditions in the area studied (Marginality < 0.9) (Hirzel et al., 2002).L. achine is recorded in areas where the climatic conditions differ greatly from the average conditions in the study area (Marginalities between 0.9 and 1.5), a fact that probably accounts for the rarity of this species.Specialization values are unique for each species, with the range of climatic conditions from one to several times lower than that in the Iberian region.However, the high values of Specialization for L. achine, P. golgus and P. nausithous indicate high constraints in terms of climatic conditions for these three species.

potential distributions
Most of the species have wide potential ranges, except L. achine, P. golgus and P. mnemosyne, for which the known distributions are restricted, and to a lesser extent E.catax and P. nausithous (Figs 1-2).In most cases the potential ranges extend into areas that are close to the where they are known to occur; occasionally conditions at more remote locations are also favourable for the species.These are locations where the species could potentially be found and should be sampled in the future.This is the case for A. isabellae, E. catax and P. nausithous.The predicted ranges for all these species are mainly concentrated on the northern fringe of the Iberian Peninsula (the Pyrenees, Cantabrian Mountains and the Iberian System, Fig. 3B).In the Pyrenees ten UTM squares have conditions favourable for up to eight of the ten species listed in the Habitats Directive: 30TXM99, 30TXN42, 30TXN52, 30TXN62, 30TXN72, 30TYM09, 30TYM18, 30TYM38, 31TBH50, 31TBH61.On the other hand, the conditions in 2,977 squares are unfavourable for these ten species (Figs 3B, 4B).

Overlap with protected Natural Areas (pNA)
All the species studied showed at least one cell matching those of the PNA (39% of the presence data are within the PNA, Figs 3A, 4A).The predictive models indicate that 866 UTM squares of 10 × 10 km (71.6%) of the PNA are favourable for at least one of the species listed in the HD.In fact, two squares potentially favourable to eight of the species listed in the Directive are mainly included in the PNA (Figs 3B and 4C).

DISCUSSION
Among the Lepidoptera listed in the annexes of the Habitats Directive, the Iberian species have ranges mainly in the northern half of the Peninsula (namely the Cantabrian Mountains and Pyrenees).This last mountain chain hosts six of the ten species studied.Interestingly, the south eastern coast of the Peninsula (Murcian-Almerian phytogeographic region, Rivas-Martínez et al., 2002) hosts the lowest number of species listed in the Directive.This is consistent with other studies based on the total richness of Iberian butterflies (García-Barros & Munguira, 1999;García-Barros et al., 2000;Romo et al., 2007).
The predictions for these species generally indicate that the areas in the surroundings of their current known distribution are also favourable for them.Thus, for the species with more restricted distributions, such as L. achine, P. golgus and P. mnemosyne, only areas in the vicinity of their known present range are likely to be favourable.This may be due to their very special ecological requirements.Changes in the natural habitats of these insects resulting from human activities may further restrict their geographical ranges and put them at risk of extinction (Freese et al., 2005;Konvicka et al., 2008).In fact, L. achine and P. gol gus are already included in the Red Data Book of Threatened Invertebrates of Spain within the categories Vulnerable and Endangered, respectively (Verdú et al., 2011).The fact that the potential distributions add a very small area to the known ranges is another cause for concern, because the possibility for broadening the ranges of these species is limited by climatic conditions.
A problem of this approach (modelling species potential areas based on climatic variables) is that the prediction of Fig. 4. Frequency distribution of the number of species on the Iberian Peninsula listed in the Habitats Directive that occur in Protected Natural Areas (A), the number of squares predicted as favourable for at least one species (B) and the number of those squares included in PNA (C).Note that the scale of the X axis is not linear but equal to the square root of the area.favourable areas is not based on habitat preferences.As a result, the areas selected by this method could no longer be suitable for the species listed in the HD.Thus, further field evidence of the species presence and their abundance would be needed for a more efficient conservation strategy.
The wide potential distributions of the remaining species indicate either areas that have or had favourable conditions for the survival of the species.However, there are nonclimatic factors that together with isolation barriers might determine the presence of butterfly and moth species, such as inter-specific competition or the lack of host plants (Schweiger et al., 2012;Romo et al., in press).Despite the aforementioned preponderance of favourable areas in the northern region of the study area, for some species the prediction may be more worrying.For instance, it is predicted that there are remarkably large areas favourable for Actias isabellae in the western half of the Peninsula where climatic conditions are suitable.However, the larvae of A. isabellae feed on pine needles (namely Scots pine -Pinus sylvestris, and Austrian pine -P.nigra).These trees are generally assumed not to be native in the NW of the Peninsula (Blanco et al., 1998), so this may account for the absence of this moth in these regions.Our prediction is also in accordance with that of other authors (Chefaoui & Lobo, 2007), who account for the current distribution of A. isabellae in terms of historical factors.Proserpinus proserpina and Eriogaster catax are under-sampled species, and therefore prediction of the areas favourable for them is exaggerated when compared to other species and their known distribution.P. proserpina flies at night in seasonally flooded riparian areas on the middle and lower reaches of rivers, a habitat greatly altered in the last few decades by human activity.The presence of favourable areas for Phengaris nausithous in Huesca and Teruel (in the southwest of the Pyrenees) also requires some comment.The persistence of this butterfly depends on the presence of its food plant (Sanguisorba officinalis) and a particular species of ant (Myrmica rubra or Myrmica scabrinodis), which it needs in order to complete its life cycle, and traditional agricultural activities that promote the development of grassland areas (Van Swaay et al., 2012).In addition, the dispersal ability of this butterfly is moderate, which makes it unlikely that individuals from the closest populations could reach favourable areas in Pyrenees (in fact an extensive exploration of the Pyrenees has failed to find this species).
All species listed in the Habitats Directive have populations within PNA.This is because although these protected areas were not created for Lepidoptera, they include the areas with greatest number of species of butterflies (Romo et al., 2007).Because of the distribution of the species it is unlikely that one cell will host all the species.For that reason, we recommend the inclusion in PNA of at least those cells with conditions favourable for eight species listed in the Directive.Of the ten squares that fulfill this criterion, three (30TXM99, 30TYM09, 31TBH50) are not currently included in NPA and should be protected.Only two squares are mainly included in the current NPA (30TYM18 and 30TYM38), while the percentage of surface coincident with a protected area for the remaining five squares is less than 15% and therefore we suggest that the protected area should be increased to completely include these five squares.

CONCLUSIONS
The squares with high densities of species listed in the Habitats Directive are located in the Pyrenees, where six of the ten species studied occur, and also areas with conditions favourable for these priority species are mainly concentrated in the northern half of the Iberian Peninsula (the Pyrenees, Cantabrian Mountains and Iberian System).Not all these areas are in the Iberian Protected Natural Areas network, and for this reason we propose the inclusion of at least eight squares with favourable conditions for most of these priority species.This study indicates that modelling the distributions of endangered species can provide an objective way of identifying potentially favourable areas for conservation and so increase the chances of survival of these species.
ACKNOWLEDGEMENTS.The study was partly funded by a "Species of Community Interest" project funded by the Spanish Ministerio de Medio Ambiente and TRASEGA.Funds were also provided by project POII11-0277-5747/BANDENCO (Junta de Comunidades de Castilla-La Mancha).We are grateful to J. Ylla for his help in determining the distribution of the moths studied.We are indebted to both reviewers whose comments substantially improved the first version of the manuscript.

Fig. 1 .
Fig. 1.Actual and potential distributions on the Iberian Peninsula of species listed in the Habitats Directive.White squares indicate current known distributions of the species.Areas predicted to have the most favourable conditions for these species are shown in black and those less favourable in different shades of grey.A -Euphydrias aurinia; B -Actias isabellae; C -Parnassius apollo; D -Phen garis arion; E -Proserpinus proserpina.

Fig. 2 .
Fig. 2. Actual and potential distributions on the Iberian Peninsula of species listed in the Habitats Directive.White squares indicate current known distributions of the species.The predicted favourableness of the different areas for each of the species is indicated on a scale black = highest and light grey = lowest.A -Eriogaster catax; B -Lopinga achine; C -Parnassius mnemosyne; D -Phengaris nausithous; E -Polyommatus golgus.

Fig. 3 .
Fig. 3.A -Number of species listed in the Habitats Directive (HD) recorded on the Iberian Peninsula.Warmer colours indicate squares with the highest numbers of these species.B -Squares on the Iberian Peninsula with conditions predicted to be favourable for the species listed in the Habitats Directive.Warmer colours indicate areas favourable for several of these species: red favourable for eight and blue for one of the ten species listed in the HD.In both A and B the Protected Natural Areas are indicated by grey hatching.

Table 1 .
Bioclimatic variables (and their abbreviations) used to model the distributions on the Iberian Peninsula of those species listed in the Habitats Directive.

Table 2 .
Bioclimatic variables that determine the distributions on the Iberian Peninsula of those Lepidoptera listed in the Habitats Directive, showing the percentage of the variation accounted for by each of the variables selected and the Marginality (M) and Specialization (S) values for each species.