Life History Traits Reflect Changes in Mediterranean Butterfly Communities Due to Forest Encroachment

The biodiversity of the Southern Balkans, part of the Mediterranean global biodiversity hot-spot, is threatened by land use intensification and abandonment, the latter causing forest encroachment of formerly open habitats. We investigated the impact of forest encroachment on butterfly species richness, community species composition and the representation of life history traits by repeated seasonal visits of 150 one-hectare sites in five separate regions in three countries—Greece, Bulgaria, and the Republic of Macedonia (FYROM—the Former Yugoslav Republic of Macedonia)— 10 replicates for each habitat type of grasslands, open formations and scrub forest within each region. Grasslands and open formations sites hosted in average more species and more red-listed species than scrub forest, while no pattern was found for numbers of Mediterranean species. As shown by ordination analyses, each of the three habitat types hosted distinct butterfly communities, with Mediterranean species inclining either towards grasslands or open formations. Analysing the representation of life history traits revealed that successional development from grasslands and open formations towards scrub forest shifts the community composition towards species overwintering in earlier stages, having fewer generations per year, and inhabiting large European or Eurosiberian (e.g. northern) ranges; it decreases the representation of Mediterranean endemics. The loss of grasslands and semi-open formations due to forest encroachment thus threatens exactly the species that should be the focus of conservation attention in the Mediterranean region, and innovative conservation actions to prevent ongoing forest encroachment are badly needed.


Introduction
The Mediterranean region of Europe is a global biodiversity hot-spot, due to its exceptional endemism rate, species richness and threat degree [1]. The high diversity of mountain ranges,

Ethics statement
The study was carried out in accordance with the national laws and permits obtained from authorised institutions: Bulgaria (National Museum of Natural History, Sofia), Greece (Eλληνική Δημοκρατία, Yπoυργείo PεριβάλλοντοB, EνέργειαB & KλιματικήB AλλαγήB, Eιδική Γραμματεία Δασών, Γενική Διεύθυνση AνάπτυξηB & PρoστασίαB Δασών & Fυσικών Pόρων, Δ/νση Aισθ. Δασών, Δρυμών και ΘήραB Τμήμα Γ & B, No. 170916/1344), and Macedonia (Bird Study and Protection Society of Macedonia). The fieldwork was not carried out in any privately owned nor protected areas. All the butterflies were carefully handled and released after identification; we collected up to five individuals per visit only for taxa not identifiable in the field, for genital preparation and species identification in the lab (not applicable regarding protected species).

Site selection
Our study area was located in the Southern Balkans, encompassing five regions (R1-R5): three in Greece, one in Macedonia, and one in Bulgaria (Fig 1), (S1 Text). We predefined three forest encroachment categories, in terms of woody vegetation cover (> 1.5m), with the help of post-2010 aerial photographs: (a) Grasslands-herbaceous vegetation dominance and a woody plant cover less than 5% with tracks of active grazing, (b) Open formations-near-even representation of woody and herbaceous cover, (c) Scrub forest, with dominance of woody plant vegetation above 70%. We located 30 sampling sites of 1ha standard area in each region (altitude from 10 to 1100 m a. s. l. (mean 440 ± 254 SE)), so as to equally represent the three forest encroachment categories, resulting in an overall number of 150 sites sampled.

Life history traits
We considered 15 life history traits for butterflies readily available in literature and reflecting (a) specialism vs. generalism (Feeding index, Flight period, Generation numbers, Migration, Overwintering stage, Wingspan), (b) larval feeding habits (Gregariousness, Host plant form, Larval feeding mode, Myrmecophily, Ovum placement), as well as (c) distribution profile (Altitudinal range, Mountain distribution, Range size, Range type). These functional traits are linked with the resilience of the species to environmental or land use change and hence its inherent vulnerability tendency (S2 Table). Information on most of the traits is directly available in literature, or, as in the case of Feeding index, easily calculable from published data (see S3 Table). An exception was information on range size, where we used simple numeric coding based on published distribution maps (see S3 Table for references).

Environmental parameters
We collected the following 15 environmental parameters for each site sampled. Two variables were collected to describe the forest encroachment gradient and were inserted in the model as predictors, namely: Forest encroachment was categorically variable with tree levels (Grassland, Open formations, Scrub forest); Canopy cover (percentage cover of woody species >1.5m) was a continuous predictor.
For each site visit, we recorded further parameters describing momentary weather conditions, namely Air temperature, Cloudiness (1: clear to 3: half-sunny), Wind (1: calm to 3: moderate breeze), as well as momentary Nectar supply, using a simple ranked scale (1: none or a few isolated flower heads, 2: isolated flowering patches, 3: whole site in bloom).

Data analysis
We transformed the recorded semi-quantitative butterfly abundances to mean numbers of individuals within the respective quantitative intervals, summed this across the four visits, and log-transformed. Air temperature, Cloudiness, Wind and Nectar were also summed across the four visits to obtain more detailed scales.
We applied generalised linear models in R [30] (Poisson distribution of the response) to analyze Forest encroachment and Canopy effects on species richness, numbers of Red-listed species, and numbers of Mediterranean species. For all three response variables, we first tested independent effects of the two primary predictors, considering also polynomial response for Canopy. Next, we tested independent effects of all site parameters and visit parameters (i.e., potential nuisance covariables), and used stepwise selection based on all potential covariables, evaluating alternative models' fits according to the Akaike information criterion (AIC) to obtain covariate models, defined as models best explaining the response variables without referring to the predictor(s) of interest. Finally, we manually forced the predictors Forest encroachment and Canopy onto the covariate models, thus assessing their effects while statistically controlling for variation due to nuisance variables.
To study changes effects on species composition, we used redundancy analysis (RDA), a constrained linear ordination, using CANOCO 5 [31]. We first computed single-term ordinations for both predictors of interest and all covariables. Next, we defined a covariate model, based on forward selection from potential covariables. Finally, we computed partial RDAs with predictors Forest encroachment and Canopy, controlled for effects of covariate model terms. We log (x + 1) transformed and centred species abundances in all RDA analyses, and evaluated significances of the ordinations using the Monte Carlo test (999 permutations).
We used the partial RDAs to analyse the life history traits responses. Because life histories co-vary with phylogeny (e.g. [32]), we constructed a phylogenetic tree of all recorded species, based on published phylogenies, supplemented by formal classification into genera and subgenera (S2 Text). We turned this tree into a patristic distance matrix, representing the distance of any pair of taxa measured along the branches of the phylogenetic tree. We transformed this distance matrix into a set of descriptors using principal coordinate analysis (PCoA), with PCoA scores centred and standardised. Not all PCoA scores are related to response variable, therefore we used their subset selected by forward selection-only descriptors with p < 0.04 were included. Finally, we interpreted the species traits responses to Forest encroachment and Canopy individually for each trait, after removing the variation explained by phylogenetic descriptors. We evaluated each step using the Monte Carlo test (999 permutations).
Raw dataset for all analyses is available as supplementary material S1 Raw Data.
When tested individually against all predictors and covariables, species richness responded to Forest encroachment, being highest in Grasslands and lowest in Scrub forest, and decreased linearly with Canopy cover (Fig 2A and Table 1). Of all potential covariables, Nectar had the strongest separate (positive) effect. Regarding site characteristics, richness was highest in intermediate longitudes, and increased with altitude. It also increased with presence of Water and Herdsman's hut. The strongest site characteristic effect, Veg1, pointed to richness increasing with humidity. The combined covariate model (Table 1) explained over 38% of variation in per site species richness. Adding the predictors of interest to this model did not reveal differences among the three stages of Forest encroachment, but revealed a significant decline with increasing Canopy (Fig 2B).
Numbers of Mediterranean species did not differ between Forest encroachment categories nor responded to Canopy. They responded to geography covariates and increased at steep slopes sites affected by grazing (the combined covariate model explained 39.8% of variation). Red-listed species, in contrast, responded significantly to both Forest encroachment (much lower in Scrub forest: Fig 2C) and Canopy ( Fig 2D, polynomial decrease). On the other hand, they did not respond to any covariates except for Nectar (3.1% of variation). After controlling for nectar, the effects of Forest encroachment and Canopy remained significant.

Species composition
In the single-term ordinations, both Forest encroachment and Canopy significantly affected species community composition. The explained variations were rather low, however, if compared with covariate predictors such as Veg1, Nectar or Longitude ( Table 2). The forward selection procedure selected the following covariate model: Altitude + Cloudiness + Latitude + Longitude + Nectar + Slope + Veg1-4 + Water + Altitude × Latitude × Longitude (36% of variation, F = 6.3, P = 0.001). On residuals of this model, both Forest encroachment and Canopy retained their significant effects ( Table 2).
In the partial RDA with Canopy (1.85%), practically all Mediterranean species, as well as practically all Red-listed ones, inclined towards low Canopy (Fig 3B).

Species traits
Visualization of relationships among the life history traits showed a clear difference between large and mobile species with multiple generations per year, overwintering in later stages,  Table 1  occurring as adults late in season, having wide host plant spectra, and inhabiting large Holarctic or Eurosiberian ranges, and species with opposite traits, typically with restricted Mediterranean or European ranges (S2 Fig). The second gradient distinguished species with multiple generations, developing on forbs and/or consuming generative plant parts, from those forming few generations per year, and feeding on woody plants or grasses.
Only a few traits responded significantly to both Forest encroachment and Canopy (Table 3). Generation numbers responded to both predictors, indicating that closed forests were inhabited by butterflies forming fewer generations per year.
Regarding Forest encroachment (Fig 4A-4D), species overwintering in earlier stages displayed affinity towards Scrub forest. Vegetation closure decreased the representation of Mediterranean and Holarctic species and increased that of European and Eurosiberian species. Spring and autumn species prevailed on Grasslands, summer species inclined towards Open formations or Scrub forest. Regarding Canopy (Fig 4E-4H), there were marginally significant

Discussion
The large-scale comparison of South Balkan butterfly communities indicated that compared with grassland and open formations, sites overgrown by scrub forest hosted lower species richness and lower richness of Red-listed species, but the same number of Mediterranean species.
In ordination analyses, we found profound changes in the butterfly community composition

Species richness and community composition along forest encroachment gradient
Forest encroachment, expressed either as a categorical predictor or as a proportion of woody Canopy cover, was associated with local butterfly richness decline regardless of the site characteristics and visit circumstances covariables for the canopy cover case. It is well known that a majority of European butterflies avoid closed-canopy habitats [33], and hence it is hardly surprising that canopy closure represents a direct threat to this insect group. Our results thus corroborate, over a relatively large geographic scale, the dependency of many butterflies occurring in the Mediterranean region on open formations (grasslands, open forests), previously documented for Mediterranean species in local-scale studies (e.g. [13,16,34,35]). Covariables increasing butterfly species richness included Water presence and the vegetation gradient Veg1, both revealing that lack of humidity restricts local species richness in the Mediterranean [36]; and Herdsman's hut, suggesting positive effects of grazing-associated disturbance on species richness. Village proximity affected species richness negatively, indicating  Table 2). See Table 2 for results of statistical tests. Species with Mediterranean ranges written in bold, Red-listed species in CAPITALS. doi:10.1371/journal.pone.0152026.g003 Forest Encroachment and Mediterranean Butterflies that species richness was not supported by other human activities than grazing. The negative effect of the vegetation covariable Veg3 (distinguishing natural and weedy communities) supported the latter conjecture. Notably, species richness increased with altitude, which seems to contradict well known patterns of altitudinal richness decline [37,38] but this was due to the fact that our sampling was restricted to lower elevations, not covering high mountains, while the biodiversity of the elevations in the Mediterranean seems to be drought restricted [36].
The richness patterns were strikingly different if only Mediterranean species or only Redlisted species were considered. For the former, we failed to detect a dependency on any of the two predictors describing forest encroachment. We also found meaningful responses of this group of species to potential covariates, although sometimes contrasting to those for total species richness (e.g., Mediterranean species increased, rather than decreased, with Temperature, and responded oppositely to the major vegetation gradient Veg1). For the latter, increasing canopy cover was by far the best predictor restricting their numbers, and the only significant covariate was (rather trivially) Nectar. These contrasting results arguably reflect definitions of the two groups. The Red-listed group contains species of all possible distribution ranges, but sharing a high degree of threat within Europe, and loss of open habitats threatens European butterflies in general [28,39]. In contrast, the Mediterranean group is defined by shared distribution range, independently of habitats, and our samples included species of all possible habitats, from bare grounds (e.g., Carchardodus orientalis, Chilades trochylus, Pseudochazara anthelea) to closed forest (Kirinia roxelana, Zerynthia cerisy) (cf. [40]). Thus, apart from the Table 3. Results of life history traits analysis. Traits-based interpretation of partial RDA ordinations of Southern Balkans butterfly community species composition (BSC) that assessed the response to Forest encroachment and Canopy models including significant covariates and controlled for phylogeny.

BSC~Response + [Trait] | Covariates
Forest encroachment Canopy Covariate model structure as in Table 2. AEV = Adjusted explanatory variable (%) Significance as follows: ** < 0.01; * < 0.05; . < 0.1.   Table 2 for formulation of covariate model) and removing the effects of phylogeny. The arrows in panels (A-D) stand for horizontal ("Scrub forest") and vertical ("Open formations") ordination axes in Fig 3(A), whereas panels (E-H) refer to ordination diagrams in Fig 3(B). Statistical tests in Table 3. doi:10.1371/journal.pone.0152026.g004 Forest Encroachment and Mediterranean Butterflies low number of Red-listed species in closed canopy sites, analysing mere species numbers does not convey much information regarding individual species requirements. The ordination analyses focusing on species composition provided deeper insights. Treating forest encroachment levels as a categorical predictor showed that each of the three categories hosted some Mediterranean and some Red-listed species, although both Grasslands and Open formations hosted apparently more such species than Scrub forest. Moreover, each of the two categories attracted distinct species, suggesting that to sustain the whole butterfly diversity associated with traditional Mediterranean landscapes, mosaics of alteration of grasslands open "savannah-like" formations are necessary. Note that even open habitat species may temporarily utilise cooler microclimates provided by close canopy sites [41], which explains our scrub forest records of such species as Hipparchia statilinus, a Mediterranean species requiring near-bare ground for larval development (27 forest presence records out of 85) (cf. [42]). Open formations hosted both Mediterranean species, some of them threatened (e.g., Hipparchia syriaca), in combination with species that prefer barren surfaces in more northerly parts of their ranges (e.g., Hipparchia statilinus, Pseudopilotes vicrama) [28,43].
The patterns found for the numeric predictor Canopy were even clearer, revealing avoidance of Mediterranean species, and affinity of northern species, towards increasing Canopy cover.

Species traits changes along forest encroachment
The prevailingly European species colonising Scrub forest sites form few generations per year, overwinter in early stages, and, counter-intuitively, fly in spring; the prevailingly Eurosiberian species colonising Open formations fly mainly in high summer; and the prevailingly Mediterranean species colonising Grasslands were mainly spring or autumn flying species. It follows that a link exists between species ranges, habitat successional stage and associated butterflies' development. This was previously suggested by Dennis et al. [21], who related the life histories of British butterflies to the life history strategies of their host plants. Association of slowly developing species forming few generations per annum with late successional habitats has been reported from such disparate regions as Germany [44], Catalonia [45] and Japan [46]. This is sometimes attributed to habitat disturbance dynamics, in that rarely disturbed habitats allow for slower insect reproductive rates in contrast to frequently disturbed habitats. This interpretation, however, fails to explain why woodland species both overwinter in early stages and occur as adults early in spring, which forces them to develop rapidly. An alternative explanation, suggested, e.g. by Cizek et al. [32] invokes the nature of antiherbivore defenses in late-successional plants (trees, coarse grasses). In such plants, quantitative defenses (tannin, silica et.) prevail, restricting associated herbivores' development to young plant tissues, available in early season. In parallel, woodland species developing on forbs are constrained to early development by rapid canopy shading, or progressive host plant senescence (e.g. [47]). The marginally significant effect of host plant form in our analysis circumstantially supports the plant defenses role. Moreover, species with higher generation numbers and species overwintering in later stages inclined towards grasslands, where the combined effects of host plants senescence and canopy shading do not apply.
Similar logic may explain the link between Mediterranean distribution, spring plus autumn adult period, and Grasslands. Grasslands receive enough sun early in the year, get hot and dry during high summer, but become inhabitable again with autumn rains [48]. Then, multivoltine species (e.g. Gegenes pumilio, Pieris krueperi, Chilodes trochylus) form additional generation(s), whereas univoltine species with long-living adults (cf. [49]) locate both nectar and oviposition substrates there. The association of species flying in high summer with Open formations is best explained by the structural heterogeneity of such sites, where mosaics of closed and open vegetation offer varying microclimate conditions, supplying some nectar, moisture and shade even during summer.

Conservation Implications
Species' ranges result from phylogenetic history, dispersal and habitat requirements [50]. The avoidance of closed canopy sites by the range-restricted Mediterranean species, and their affinity for either Open formations or Grasslands, agrees with results recently reported for Greek birds [51] and spiders [14]. Grill et al. [16] and Kati et al. [12] reported, for butterflies and orthopterans, respectively, the highest species richness, and highest representations of rangerestricted species, from such richly structured habitats as abandoned orchards and wooded pastures in the Greek nature reserve Dadia. Increases of common northern species at the expense of Mediterranean endemics were also detected for southern French birds [52], Sardinian plants [53], and Catalonian (i.e., West Mediterranean) butterflies [15]. Assuming historical conservatism of species life histories, the negative association of Mediterranean species with closed canopy condition falsifies the "forested Mediterranean" hypothesis, highlighting the need to maintain open landscapes across the region.
Notably, the increases of northern species due to forest encroachment contradict the predictions that northern species should decline at their southern range margins due to the current climatic warming [54]. This process is probably counteracted by another development, detected for Greek butterflies by Zografou et al. [55], who found an increase of low-altitude thermophilous species against high-altitude ones. The two processes, increase in the representation of thermophilous species due to warming climate and their decrease due to habitat loss, are likely affecting species individually, depending on their ability to adapt, e.g. by locating sites with suitable microclimates [56]. For global conservation, however, the outcomes are hardly positive, because the majority of the Mediterranean endemics depend on grasslands or open formations, habitats that are rapidly decreasing all over the study region.
Without maintaining rich mosaics of open and semi-open habitats across the southern Balkans, the restructuring of butterfly and other small animal communities due to forest encroachment will gradually replace range-restricted endemic fauna by wide ranging generalists. Maintaining open landscapes is complicated by several factors. First, such widely advocated land management tools as "headage payments" for shepherds [8] or agro-environmental schemes rewarding environmentally benign farming [57], were originally designed in northwestern Europe and may be poorly transferable to the conditions of Southern Europe, with much more diverse habitat conditions and declining rural population [58]. Second, financial incentives do not guarantee that human impacts on habitats replicate those existing in the past. For instance, agrotechnology developments such as fodder crops production and vehicle transport relaxed the need to harvest summer coppice, or to move herds across the landscapes (transhumance) (cf. [59]). Cizek et al. [60] documented that current management technologies fail to provide microhabitat heterogeneity needed for reserve management in Central Europe, and the outcomes may be even worse in species-richer Southern Europe. Still worse, relying on subsidies assumes constant economic growth, which is far from guaranteed in the long term. Economic decline might promote returns of urban population to villages, but this would be a long-term process, whereas breakdowns in funding may lead to rapid habitat and species losses.
Without downplaying the subsidised efforts to maintain rural habitats diversity [61], novel approaches which would maintain the open to semi-open conditions across the Mediterranean while being economically sustainable should be sought. At least locally, declining grazing by farm animals might be replaced by free ranging ungulates, including species historically extirpated from the Mediterranean [62]. Such projects would, at least regionally, return to the Mediterranean the key players that had been affecting ecosystem dynamics before the advent of farming, and with which the regionally endemic biodiversity has evolved.