Environmental factors influencing butterfly abundance after a severe wildfire in Mediterranean vegetation

Environmental factors influencing butterfly abundance after a severe wildfire in Mediterranean vegetation.— Despite the attention given to the ecology of butterflies, little is known about their community response to wildfires in the Mediterranean region. Here, we evaluated the butterfly assemblage two years after a severe, 13,000 ha wildfire in Catalonia (NE Spain) in relation to the surrounding unburned habitat. Using visual transect censuses we assessed community parameters such as abundance, diversity, species richness and equitability in burned and unburned areas. Correspondence analysis was used to analyse specific composition and relative abundance of species in the community. The influence of environmental variables on the abundance of some common species was analysed using generalized linear mixed models, taking spatial effects into account. No significant differences were found between areas for any of the community parameters, and no dominance was detected in the burned area. The structure of the vegetation and the geographical distribution of transects influenced the ordination of species and transects on the correspondence analysis plot. Generalized linear mixed models (GLMM) results underscored the role of nectar availability, fire and vegetation structure on the abundance of most species studied.


Introduction
Wildfires are a major ecological disturbance, affecting ecosystem functioning and species composition in forests around the world (Bond et al., 2005;Blondel et al., 2010;Mateos et al., 2011). In the Mediterranean region in particular they are considered an ecological factor that shapes the landscape and the ecosystems (Lloret, 1996). The risk, frequency and intensity of forest fires have increased, however, in recent decades largely due to land use changes. Pastures and agricultural land, for example, have been abandoned in mountain areas, and tree plantations have been established for commercial use. Furthermore, the cessation of traditional forestry has led to fuel accumulation over large areas (Debussche et al., 1999;Feranec et al., 2010) Climate change has also contributed to an increase in the frequency of fires, and is expected to cause even greater impact in the future (Piñol et al., 1998).
In a context of global change, studying the response of species to environmental disturbances has become crucial and one of the main goals for conservation and landscape management (Bengsston et al., 2000). More particularly, in the Mediterranean and other regions, the responses of invertebrates to fire have been examined in diverse taxonomical groups as a way to quantify the effects of fire on species distribution and abundance, and also on changes occurring at the community level (Swengel, 2001;Kiss & Magnin, 2003;Moretti et al., 2004;Santos et al., 2009;Mateos et al., 2011). However, quite surprisingly, little is known about how butterfly communities are affected by wildfires, even though this group is considered an excellent indicator of biodiversity trends in terrestrial ecosystems (Thomas et al., 2004).
Butterflies have highly specific requirements in terms of feeding resources in both the larval and adult stages (Erhardt & Mevi-Schutz, 2009;Munguira et al., 2009) and regarding the microclimatic conditions needed for the viability of populations (Thomas et al., 1999;Roy & Thomas, 2003). Some species have limited mobility and live in meta-populations, being strongly and negatively affected by habitat destruction and landscape fragmentation (Hanski & Thomas, 1994;Steffan-Dewenter & Tscharntke, 2000;Bergman et al., 2004). All these features and the ease with which they can be monitored make them an ideal target to explore the effects of forest fires on terrestrial insects.
In this paper we report a study carried out to document the response of butterfly communities in an area in NE Spain that was severely affected by a forest fire. The analysis is presented at two levels: first, at the community level, to describe the effects of fire on the composition of the communities studied, and second, at the species level, to investigate the main factors affecting the relative abundance of some common species within these communities.
A large wildfire like the one we examined is a first-order disturbance on the flora and fauna as it drastically reduces food resources and causes massive mortalities of organisms. It can also have indirect effects on the structure and species composition of plant and herbivore communities. It has been reported that various insect species (including some butterflies) decrease sharply in numbers in the early stages after a fire (Swengel, 2001). Because of their sensitivity to environmental alterations and changes in vegetation structure, butterfly communities are strongly affected by forest fire. We predicted a decrease in butterfly diversity in burned areas due to the local extinction of some species and a greater dominance by opportunistic species able to recolonize the primary successional stages after such a severe disturbance (Odum, 1969;Steffan-Dewenter & Tscharntke, 1997).
The abundance of a butterfly species is determined by a combination of environmental and biological factors. We expected to find a strong influence of factors related to the availability of food resources (nectar availability and abundance of larval host plants), vegetation structure (i.e., cover of different vegetation layers) and the fire effect itself. The importance of these factors depends on specific biological characteristics that determine the sensitivity of species and the resilience of populations.
In connection with these predictions, the goals set in this study were: (1) to assess the modification of the butterfly community two years after a large wildfire in relation to control areas, using general descriptors such as species richness, abundance and dominance; and (2) to analyze the abundance of some common species and determine the most influential environmental factors in the recovery of butterfly populations.

Study area
The study area is located in the county of Alt Empordà (Girona province, NE Spain) ( fig. 1). The region has a rugged relief and the climate is subhumid and humid Mediterranean, with the strong influence of the northern wind, known as the tramuntana. Although the potential vegetation is holm oak (Quercus ilex) forests, current landscape is a mosaic resulting from historical and current land use, with a dominance of Aleppo Pine (Pinus halepensis) and abandoned cork forests and crops (vineyards and olives) that have now turned into Mediterranean shrub land. The fire in the region occurred on 22nd July 2012 and lasted six days. Driven by the tramuntana, 13,963 ha (according to the Forest Division of the Catalan Government Fire Service 2013) were burned. The region has a long history of wildfires but this was the largest since 1986 when 26,000 ha burned. Post-fire management in the 2012 fire consisted mostly of logging, with timber being removed or made into chips for use as fuel for power plants.

Sampling design
We selected seven sampling localities. To reduce environmental variability, all localities were situated in pine forests and shrublands on the western part of the burned area and its nearby unburned area, on

Butterfly sampling
Butterflies were sampled every two weeks from the beginning of April 2014 until the end of June 2014, with a total of five visits to each transect. The sampling method consisted of counting adult butterflies following the standards of the Butterfly Monitoring Scheme (BMS) method (Pollard & Yates, 1993). Transects were walked at a constant speed, between 11am and 5pm, under appropriate weather conditions (> 50% sun, wind ≤ force 3 in the Beaufort scale, temperature ≥ 17ºC). We counted only those butterflies that were 5 m or less in front and 2.5 m to the sides of the recorder. When identification to species level was not possible at first sight, butterflies were captured with an entomological net for close inspection and then released. Butterflies were identified consulting Tolman & Lewington (2011) field guide.

Species abundance modelling
A subset of 15 butterfly species was selected to study the influence of environmental and biological factors on their relative abundance. Species were selected mainly based on their commonness at the study sites, and also because their flight period coincided with our sampling period. Table 1 gives a list of the 15 species with some basic information on their natural history in the study area and the response of their host plants to fire disturbance. We also measured the availability of food resources (nectar for adults and host plants for larvae) and vegetation structure (foliage cover of the different strata) to model butterfly abundance in each locality (cf. Dennis, 2010).
Foliage cover was estimated once for three vegetation strata (grass, shrub and tree) comparing it to a standard template (Prodon & Lebreton, 1981). For each transect, foliage cover was estimated at three equidistant sites and the mean was calculated.
Nectar availability was estimated in each census using a semi-quantitative scale of flower abundance.   For the first four plant species, abundance was directly measured by counting the number of individual plants along the transect. This was not possible for the other plant species due to their high density. Abundance was then estimated from their cover using the semiquantitative scale: 0. Species absent (0%); 1. Low cover (< 25%); 2. Moderate cover (25-50%) and 3. High cover (> 50%).

Data analysis
Butterfly data were summarized in a specific composition table showing total abundance of each species for burned and unburned control transects. Data were analysed using the R statistical software (R CoreTeam, 2014). First, we conducted a comparative analysis of the structure of the community in burned and unburned areas. A correspondence analysis (CA) implemented with the statistical package BiodiversityR (Kindt & Coe, 2005) was used to assess the structure of the community and to explore the main gradients influencing butterfly composition. The input for the CA was the total count of each species per transect. In this analysis, we excluded species whose total count was less than five individuals.
Second, to highlight differences between the two treatments (fire vs. control), we also calculated the following community descriptors: richness (S), abundance, diversity (Shannon-Wiener index, H) and equitativity (Evenness Index, E): Normality and homoscedasticity of the descriptors was tested using Shapiro-Wilks and Bartlett tests. Differences in the descriptors between the two treatments were tested using ANOVA.
The abundance of the most common species was modelled using generalized linear mixed models (GLMM) to assess the importance of four environmental variables: fire effect, vegetation structure, and availability of food resources for adults, and for larvae. Six potential explanatory variables were selected as descriptors of these environmental variables and were used as fixed factors on the GLMM: fire (Fire), percentage of foliage cover of the three vegetation layers: herbaceous (Herb), shrub (Shrub) and tree (Canop), host plant availability (Host Pl) and nectar availability (Nectr). We used locality as a random factor to account for our particular sampling design (with burned and control transects paired in specific localities) while controlling for possible site-based differences. Box-plots were used to check the potential influence of fire on the environmental variables. Variables highly influenced by fire were excluded to avoid redundancy (CANOP and some host plants: Quercus coccifera, Dorycnium pentaphyllum, Cistus sp.). To obtain the model that could best explain the abundance of selected species (response variable) we performed multiple GLMM for each species. Competing models were compared using Aikake's Information Criteria (AIC) based on maximum likelihood. The model with the lowest AIC value was the approximation that best fitted the data. Differences (Ai) between the AIC value of the best model and the AIC value for each other model were used to assess model performance. Models with Ai values lower than two are considered to be essentially as good as the best approximating model (Symonds & Moussalli, 2011). Analyses were performed with the statistical package lme4 (Bates et al., 2014), using the loglink function and structure of negative binomial residues.

Community level
A total of 918 butterflies belonging to 47 species were observed in the censuses. We observed 398 butterflies belonging to 39 species in burned transects, and 520 individuals belonging to 37 species in control transects (table 2). Regarding the similarity of species between treatments, 28 were common to the burned and control transects, eight were found only in the control transects (e.g., Pararge aegeria and Melanargia lachesis) and 10 were found only in the burned transects (e.g., Vanessa cardui and Coenonympha dorus). The most abundant species in the control transects was Pyronia bathseba, representing almost 16% of the total individuals. Other abundant butterflies were Gonepteryx cleopatra (14.6%) and Lysandra hispana (10.8%). These three species represented 41.4% of abundance in controls, and the 10 most abundant species attained a figure of 78.3%. In burned transects, G. cleopatra was the most common species (16.3%), followed by L. hispana (11.3%) and Satyrium esculi (9.3%). These three species represented 36.9% of the total number in the burned transects, and the 10 most abundant species represented 72.6%.
Correspondence analysis (CA) graphics show the ordination of the butterfly species ( fig. 2A) and the 14 transects ( fig. 2B) based on specific composition and relative abundance of species in transects. The first three axes explained 57% of the variance in the dataset.
Despite some overlap in the centre of the graph, burned and control polygons were segregated along a diagonal gradient in the biplot of the first two axes, going from Vanessa cardui -a species only found in burned transects and with extreme negative coordinates-to P. aegeria -a species only found in control transects and showing extreme positive coordinates.
This ordination can be interpreted as the effect of fire on the butterfly community structure, but it does not show a strong influence of this factor. The proximity of paired transects (corresponding to the same locality) indicates the influence of the spatial distribution. This was true for all sites except localities 4 and 6, which were characterised by dense forest and had very poor butterfly communities, strongly dominated by P. aegeria, a forest species. This particularity increased the distance of both sites from their paired burned transect in the ordination plot. The polarity of the burned transects polygon was given by another diagonal gradient (perpendicular to the first one, referred above) on which the association between paired transects and the differences related to the geographical distance are most obvious. At the bottom of the gradient we find localities 2, 3, 4 and 5 (at a close distance to each other and characterised by mixed forests with predominance of pine), and at the top we find localities 1, 6 and 7 (more isolated from the rest and characterised by mixed forests with predominance of oak). The three localities at the end of the gradient are also those close to hill tops where P. machaon, I. podalirius and M. occitanica abound. These three species typically show hill-topping behaviour, by which males congregate in topographically elevated points to where females fly for mating.

Fig. 2. Ordenación de la comunidad de mariposas (A) y transectos (B) en el diagrama de dispersión biespacial de los dos primeros ejes del análisis de correspondencias (CA). (Para las abreviaturas de las especies, véase arriba.)
CA2 control), whereas species with highly positive CA1 values (e.g., P. aegeria, Euphydryas aurinia, Gonepteryx rhamni, Limenitis reducta) were only found, or found in higher abundance, in control transects with dense forest. Species located in the plot centre do not seem to have any distribution pattern associated with the degree of opening of the habitat. Despite these differences in the composition of butterfly communities, we did not find significant differences in the general community descriptors between burned and unburned areas (table 3): ANOVA tests for abundance (F = 0.94, df = 1,12; p = 0.351), species richness (F = 0.348, df = 1,12; p = 0.566), diversity (F = 0.471, df = 1,12; p = 0.505), and evenness (F = 1.527, df = 1,12; p = 0.24). Table 4 summarizes the models that best explained the abundance of the most common species in the study area following AIC criteria (AIC values on appendix 2). Models were obtained for 14 of the 15 species initially considered (table 1). It was not possible to build a model for M. lachesis due to an error related to covariance of data as the result of the multiple zeros in this species.

Species level
Results included only a single best model for 12 species. However, in three cases, the best model was found within a set of two (Gonepteryx spp. and Gonepteryx cleopatra) or three models (Lysandra hispana). In total, 19 models were considered.
Concerning the importance of fire as the only factor influencing butterfly populations, just in two species (G. cleopatra and L. hispana) this variable was sufficient to explain the abundance of butterfly species. In the other two cases with the same result (C. rubi and E. aurinia), this was the only model that could be performed due to an error related to the large amount of zeros when adding more variables. It should be noted that another model is to be considered in L. hispana (Fire + Nectar).
In most species (10 out of 14), the best model (or one of the best models) was that which included the effect of the fire and nectar availability. In four species, the model including fire effect, nectar availability and structure of vegetation was the best model. In one species (I. podalirius), the best model also included larval host plant.
The effect of these variables differed depending on the species of butterfly. Fire had a negative influence on 15 of 19 models (12 species). Nectar availability had a positive effect in 11 of the 13 species where it was present. Related to the structure of vegetation, herbaceous cover had a positive influence on three of five models and shrub cover in four of five. Canopy cover was not included in the models as it showed a high negative correlation with fire, explaining the high negative coefficient of fire in the model for P. aegeria, a forest species.
Lastly, host plant abundance was irrelevant except for I. podalirius, for which it had an expected positive effect. It should be noted, however, that for nine out of the 15 species of butterflies, this variable was not used to construct the model because we did not have direct measures of host plant abundance, or it showed a high influence of fire.

Butterfly communities in the burned area
Our results did not show any significant difference between control and burned transects in terms of mean abundance, species richness, diversity, and evenness of the butterfly community. This finding contradicts our predictions of a decline in diversity in the severely disturbed burned areas (cf. Odum, 1969;Caswell, 1976), presumably associated with the local extinction of some species (i.e. those most vulnerable because of their low mobility or because the fire occurred when they were in the critical egg or larval stages) and an increase in the dominance by rapidly colonizing species. On the contrary, our data showed that many species were able to recolonize the burned area (or resisted in it) within two years of the fire. Moreover, in our study region, butterfly communities in burned and unburned sites were both dominated by a few species, leading to similar evenness values. Thus, at the unburned sites, P. bathseba -a sedentary species-was exceedingly abundant, while at burned sites G. cleopatra -a highly mobile species-reached comparable high numbers.
The speed at which butterfly communities can recover from a forest fire was noted by Nel (1986), Table 3. Community parameters: abundance, species richness, Shannon Index and Evenness (mean ± SD).

Abundance
Richness Shannon Index Evenness Burned 56.9 ± 42.9 16.0 ± 5.0 2.35 ± 0.33 0.69 ± 0.69 Unburned 74.3 ± 20.5 13.9 ± 8.2 1.94 ± 0.98 0.64 ± 0.64 who monitored butterfly assemblages in an area in southern France that was devastated by a wildfire of characteristics similar to those studied here. Nel (1986) recorded species occurrence during the four years after the fire and compared the changing butterfly assemblages to that known to occur in the area prior to the disturbance (Nel, 1982). The fire occurred at the end of July, and two months later the number of butterfly species was very low (six species, 10% of the initial number). However, in the following years, the recolonization process was fast, with 60% of species recorded in the second year and 80% in the third year. Nel's results coincide with ours in the Alt Empordà, as we recorded a similarity of 62% in specific composition between burned and control transects two years after the disturbance (see below). Moreover, preliminary observations from another burned area in Les Gavarres (Girona, NE Spain) suggest a similar pattern of recolonization. In this latter case, the fire took place in March 2014 and burned 359 ha; as in Nel's (1986) study, our data show an initial phase lasting a few months with no butterflies at all followed by the appearance of six species (10% of species richness in control transects) by the end of the first summer. However, although the differences in community descriptors between burned and non-burned transects were non-significant due to the speed of the recolonization process, the direct or indirect effects of fire were detected in the species composition and their relative abundances. Of 47 species found, 10 were found only in burned transects, eight only in control transects, and 27 were common to both, with a similarity of 62% in terms of specific composition. These differences also became evident in the plot of the first two CA axes, which showed some segregation of burned and control areas. Some of these differences relate to large differences in habitat structure in sampling sites, as was the case between paired transects in localities 4 and 6. Controls 4 and 6 sampled dense forest, which resulted in low densities of a few forest specialists. Closed forest in the Mediterranean region typically holds low density and species richness of arthropods (e.g., Mateos et al., 2011;Verdasca et al., 2012), in contrast with more open areas that provide a high concentration of nectar sources and attract adults of most butterflies and other insects (Jubany & Rovira, 2000). Although the sampling design sought to reduce environmental variability of the study areas, habitat heterogeneity, the need to keep a short distance between transect pairs, and the availability of severely burned areas, meant that this was not always possible. The position of our transects was the result of the trade-off between proximity (to reduce environmental variability) and distance (to avoid the border effect) between paired (burned and unburned) transects, resulting in a distance of 200-700 m from the fire perimeter.
However, differences between burned and unburned sites may also be related to the functional groups in each area. For example, Kwon et al. (2013) analysed the Lepidoptera communities for five years after a fire in Korea. At first they noted a reduction in the number of specialists in the disturbed areas compared to nearby unaffected areas, but they observed that this difference disappeared by the end of the study period when the proportion of butterfly functional guilds (i.e., generalist and specialist, based on larval host plant use and adults habitat) had returned to original levels. Nel (1986) and Cleary & Genner (2004) obtained similar results. They also noted that during the process of butterfly recovery in burned areas, the first butterflies to arrive were generalist species and that these were replaced by specialist species over the following years. Swengel (1996) and Vogel et al. (2010) also found that butterfly specialists took three to four years to recover from fire disturbances in several open areas, possibly linked to the process of colonization of the area after local extinction. This was not the focus of the present study so we are unable to exclude a role of local resistance in addition to colonization after fire to explain the results found. Several recent works point to a correlation of various life history traits that allows a species to allocate along an axis from extreme generalism to extreme specialism (Carnicer et al., 2013;Dapporto & Dennis, 2013). In this context, the work by Carnicer et al. (2013) -based on data from the Catalan Butterfly Monitoring Scheme which includes several sampling sites in Alt Empordà-offers an excellent framework to investigate this issue further.

Main factors affecting the abundance of the most common species
To model the butterfly abundance of particular butterfly species in burned and unburned areas, we constructed generalized linear mixed models that took into account those environmental factors that presumably had the strongest effects on the populations. Besides the availability of feeding resources and the general habitat structure, we explicitly tested the importance of fire disturbance in explaining butterfly abundance. Models showed the outstanding importance of nectar availability, which had a significant positive influence in 11 out of 15 successfully modelled species. This result is not surprising, as many studies have shown the key role of nectar availability in explaining butterfly distribution and abundance (e.g., Loertscher et al., 1995;Schneider et al., 2003) in temperate areas. In this respect, fire may have indirect positive effects on some mobile butterfly species, as massive blooms of some flower species (e.g., Galactites tomentosa, Cistus monspeliensis, etc.) are highly characteristic in our study area one or two years after a fire disturbance (Pons & Prodon, 1996). The dominance of the highly mobile Gonepteryx spp. species at the burned sites, where its host plant Rhamnus alaternus was found in low abundance or completely absent (in agreement with Paula & Pausas, 2013), was probably explained by this fact, as population movements in search of nectar sources are a common phenomenon in our region (García- Barros et al., 2013). However, the importance of nectar availability was not only detected in well-known highly mobile species, such as Anthocharis cardamines and Colias crocea (e.g., Stefanescu, 2000;Kuussaari et al., 2014), but also in sedentary species such as Pseudophilotes panoptes where it appeared quite unexpectedly as the single determinant of butterfly abundance.
On the other hand, host plant abundance entered as an explanatory variable in only one species, Iphiclides podalirius, with the expected positive effect. For S. esculi and C. rubi, the main host plants (Quercus coccifera and Cistus spp., respectively) were not included in the model even if we had measurements of their abundance because they were strongly and positively associated with fire. These two butterfly species were the only ones showing a positive effect by fire, possibly explained by the high densities that their host plants attain in burned areas as a result of their quick resprouting (Paula & Pausas, 2013).
Interestingly, for all other species, fire showed an invariably negative effect, which was significantly detected in almost all the modelled species. This result indicates that the recovery of butterfly populations after a wildfire event may take, in many cases, more than two years. This seems to be specially the case of Pyronia bathseba, the dominant species in unburned areas, which was markedly rarer at the burned sites, and P. aegeria, a forest species which was only found at the unburned transects.
Finally, it is worth mentioning that some models included the factors related to the structure of vegeta-tion (cover of herbaceous plants and shrubs), a result consistent with the many studies pointing to the key role of habitat structure in determining butterfly preferences (see Dennis, 2010, and references therein).
To conclude, generalized linear mixed model results evidenced the influence of the availability of trophic resources and habitat structure on butterfly abundance, but also the importance of fire as a depressor of population levels in many species. However, as suggested by Cleary & Grill (2004), many other factors play direct or indirect roles in determining the presence and abundance of butterflies after a fire. These factors can be purely environmental (e.g., changes in humidity that affect sensitive species such as some Satyrines in forest habitats: Hill, 1999), or ecological (such as alterations of complex interactions between various species as in the case of myrmecophilous lycaenids). Undoubtedly, these additional factors influence the recovery of butterfly populations after a wildfire event and may account for the relatively low explanatory power of our models.