EUROPEAN JOURNAL OF ENTOMOLOGY EUROPEAN JOURNAL OF ENTOMOLOGY Ibero-Supra-Distribution of butter ﬂ ies (Lepidoptera) in a successional mosaic of Mediterranean mountain habitats

. The distribution of butter ﬂ ies was analysed in a mosaic of different types of habitats in a typical Ibero-Supramediterra-nean plant landscape. This mosaic landscape is composed of oak forests ( Quercus pyrenaica ) and their corresponding shrub and grassland successional communities. The observed patterns were based on butter ﬂ y- ﬂ ower interactions in two consecutive years (2017–2018) recorded in different habitats. The results of the nestedness analysis indicated that the same butter ﬂ y community exploits all of the successional plant-communities, but some differences due to the availability of ﬂ owers. The foraging for nectar sources was mostly restricted to a few ﬂ owering plants: Rubus ulmifolius , Carduus carpetanus , Thymus pulegioides and Dianthus deltoides . Some butter ﬂ ies changed their ﬂ ower-visiting patterns over the two years studied. The distribution of butter ﬂ ies was strongly linked to the distribution of their preferred nectar sources and changes in the use of these sources modi ﬁ ed habitat use. The area of Rubus shrubland is a small but highly used habitat in this area, given the presence of Rubus ulmifolius , a plant species with high coverage and abundance of nectar, makes this area much better for foraging than other habitats. The Cytisus and Rubus shrublands were used more than grassland, indicating that seral habitats can be more valuable than traditional open grassland and forest habitats for butter ﬂ ies.

The butterfl y-plant species interactions are always focused on one or very few nearby habitats without assessing all the habitats comprising an ecological secondary succession. However, this kind of information would be very useful for conservation and other ecological studies, especially in the Mediterranean region, due to the diversity and patchiness of its habitats (Slancarova et al., 2016) and confi rmed positive signifi cant responses of biodiversity to fragmentation (Fahrig, 2017). Secondary succession, which is infl uenced by pre-existing communities, is characterized by a series of successional communities that follow each other until reaching a more developed seral stage (Rivas-Martínez et al., 2011). Although forest is generally considered the most developed seral stage, this depends on local conditions and in extreme regions it could be scrubland or grassland. This study focused on Sierra de Guadarrama (central Spain), where over 50% of the Iberian butterfl y richness can be found (Gómez de Aizpúrua, 1987) and there is great diversity in other biological groups, such as fungi (Campoamor & Molina, 2001), vascular plants (Lobo, 2001), dung beetles (Martín-Piera et al., 1986) and vertebrates (Rey Benayas & de la Montaña, 2003). The most developed natural vegetation in this area corresponds to Pyrenaean oak (Quercus pyrenaica Willd.) Ibero-Supra-

INTRODUCTION
Butterfl y imagoes use nectariferous fl owers as food sources. The uses of nectar sources by butterfl ies in different habitats indicates they tend to prefer open habitats, such as grasslands, to closed ones, such as forests (Loertscher et al., 1995;Grill et al., 2005;Kitahara et al., 2008). Grasslands are perhaps the habitats where most butterfl y distribution studies in relation to nectar sources are conducted (Jennersten, 1984;Severns et al., 2006;Ezzeddine & Matter, 2008). A direct relationship is recorded between the occurrence of nectariferous fl owers and butterfl y richness in grasslands (Pöyry et al., 2009;Krämer et al., 2012), although other evidence does not show such a clear relationship (Sharp et al., 1974). When grasslands are encroached by shrubs or trees, there are different, negative and positive, responses in terms of butterfl y richness (Holl, 1995;Balmer & Erhardt, 2000;Söderström et al., 2001;Lien & Yuan, 2003;Öckinger et al., 2006). Specifi cally, a positive relationship between nectariferous fl owers and butterfl ies at the edge of forests (van Halder et al., 2011). However, most of these studies are for temperate Europe and are scarce in the Mediterranean biodiversity hotspot region (but see Marini et al., 2009;Slancarova et al., 2016;Ubach et al., 2020;Colom et al., 2021). bulbosa grassland and Nardus stricta meadows. Identifi cation of vegetation was based on fl oristic composition and physiognomy of the vegetation and phytosociological syntaxonomy (Molina Abril, 1994). A minimum of six plots was selected per habitat in the fi eld and subsequently confi rmed by aerial photography. The plots were located at least 100 m apart from each other to avoid overlaps. The number of plots was selected according to habitat type representativeness (frequency and extension) in the landscape up to a maximum of 11 plots per habitat. Thus, 11 plots of forests, 10 of Cytisus oromediterraneus shrublands, 6 of Rubus ulmifolius shrublands, 9 of Poa bulbosa grasslands and 8 of Nardus stricta meadows were studied. Plot size was 100 m 2 as this is the size commonly used (Chytrý & Otýpková, 2003) to reasonably represent the fl oristic variation in each habitat. Thus, the area studied constituted a total of 44 plots corresponding to 4,400 m 2 . The surface of every habitat in the territory was calculated using the program QGIS 3.10 (QGIS Development Team, 2014).
In each plot, the vegetation was recorded, visually estimating the percentage cover of each species of plant in each plot (100 m 2 ). Cover classes were not considered, but percentages were calculated as close as possible and based on several measurements. Subsequently, the observed coverages were corrected based on their representation in the area. But terfl y occurrence was recorded for 20 min per plot, during which the species that were found and moved within the plot were identifi ed. The recorder moved within the limits of the plot trying to count all the butterfl ies without disturbing them. The number of feeding interactions between butterfl ies and plants was recorded as the total feeding times for each species of butterfl y on the available fl owers, along with sightings of butterfl ies simply passing by.
Rec ordings were carried out in the fi rst fortnight of July in the years 2017 and 2018, on sunny days with no appreciable wind (Pollard, 1988). This period was selected to avoid signifi cant variations in the phenology and abundance of butterfl ies and plants (Pollard & Yates, 1993). The order of observation (day and time) of each plot was previously randomised. Recordings started at 11:00 a.m. and fi nished at 4:00 p.m. as this time slot is considered the optimal time for butterfl y activity (Pollard & Yates, 1993 mediterranean forests at intermediate altitudes (900-1,400 m a.s.l.). These forests have traditionally been removed to encourage extensive pasture lands, but the abandonment of traditional uses and grazing by livestock has favoured the increase of shrubs and the growth of forest in recent years (Sánchez-Mata et al., 2017).
Our working hypothesis is that the habitat for nectariferous plants describes butterfl y distribution at a local scale. In this research, for the fi rst time, a comparative study of butterfl y-fl ower interactions is carried out across an entire mosaic of habitats of an ecological successional series. As the habitats in the same ecological successional series intermingle at a local level, it is necessary to study the butterfl y community as a whole in order to analyse if there is any type of pattern present in their distribution. Our specifi c objectives involved answering the following questions: (a) How is the community of Lepidoptera structured across the mosaic of seral vegetation throughout the secondary succession? (b) What are the main interactions between Lepidoptera and fl owering plants? and (c) How does type of habitat infl uence butterfl y distribution? We believe that studies like this can be valuable for butterfl y conservation and faunistic studies of one of the most threatened groups of insects in Europe, especially in the Mediterranean area.

MATERIALS AND METHODS
This study was carried out in the vicinity of the town of Mirafl ores de la Sierra (North of Madrid, Central Spain). The area studied was located 3 km to the NW of the town, in the west Guadalix valley of Sierra de Guadarrama (Fig. 1). There is a small altitudinal gradient in this area (1,319-1,550 m, with most of the plots between 1,400 and 1,500 m). The soils are characterized by being acidic, with siliceous materials, such as granite, gneiss, schists and quartzites. It is a typical Ibero-Mediterranean mountain landscape, which includes patches of remaining natural oak forests, together with extended meadows and grasslands that support an important livestock population, mainly cows, as well as pine plantations (Sánchez-Mata et al., 2017). According to bio geo graphic regionalization, the area studied belongs to the Mediterranean Region, West-Iberian-Mediterranean Province, Carpetanian-Leonese Province, Guadarramean Sector (Rivas-Martínez et al., 2014). The area is located within the Supra-Mediterranean altitudinal belt. Data was recorded within this altitudinal belt to minimize climatic variability, which can infl uence the occurrence and abundance of butterfl ies and plants (Pollard, 1988). The natural vegetation in the area corresponds to oak forest of Quercus pyrenaica with seral shrubland communities, including mesic shrublands characterized by Gentisates, especially Cytisus oromediterraneus Rivas Mart. et al. and hygrophilous thorny shrublands dominated by Rosaceae species, mainly Rubus sp. pl., such as Rubus ulmifolius Schott in Isis (Oken). Seral grasslands may also be classifi ed as dry grasslands, characterized by Poa bulbosa L., and wet grasslands (meadows), characterized by Nardus stricta L. (Rivas-Martínez & Canto, 1987;Rivas-Martínez et al., 2011).

Field recording
A stratifi ed recording method was used to obtain a proportionate sample of habitat diversity. Five types of vegetation were identifi ed in the area studied: Quercus pyrenaica forest, Cytisus oromediterraneus shrubland, Rubus ulmifolius shrubland, Poa

Statistical analyses
To determine whether the effort spent recording was suffi cient and the number of species recorded approached the real number of species, a species accumulation curve was generated (Gotelli & Colwell, 2001). The curve was smoothed to avoid temporal or spatial biases and a non-parametric estimation of the potential richness based on Chao1 (Colwell, 2000) was calculated using the EstimateS program. The species accumulation curve was adjusted to the asymptotic curve of Clench using the CurveExpert 1.4 program (Hyams, 2010). Data were also analysed using nesting techniques to compare the plots on the basis of the number of species of butterfl ies recorded. These analyses were carried out using the Nestedness Temperature Calculator program (Atmar & Patterson, 1995. In addition, a cluster analysis based on the occurrences of butterfl y species per plot was done in order to check for possible community structure. It was calculated using the Bray-Curtis similarity index and the group average amalgamation rule by means of the Primer-e program (Clarke & Gorley, 2006).
Another cluster analysis was carried out based on the cover of each species of plant present in each plot. The cover values were log-transformed and a resemblance matrix was obtained based on Euclidean distances. An amalgamation strategy based on the unweighted pair-group average rule was carried out to fi nally obtain a distance cluster for each period recorded.
A Mann-Whitney U test was used to determine the differences between habitats in terms of the mean species richness of butterfl ies, in the two years considered. A Wilcoxon matched pairs test was also carried out to contrast differences between years within the same habitat. To verify whether there was a positive or negative relationship between butterfl y-fl owering plant species interactions, Generalized Linear Models were used. Only butterfl y species for which at least 10 interactions with plants were recorded were considered (data could be considered anecdotal below 10). The number of interactions of each butterfl y species per plot was considered to be the dependent variable, whereas the plant cover per plot of each species of fl owering plant selected by the butterfl ies, temperature and humidity were considered to be independent variables. In order to explore the possible curvilinear relationship (Austin, 1980) between the dependent and independent variables, the statistical signifi cance of the linear (y = a + bx), quadratic (y = a + bx + cx 2 ) and cubic functions (y = a + bx + cx 2 + dx 3 ) of the independent variables were determined. The goodness-of-fi t of the models was evaluated by means of the deviance statistic, which was compared by means of F tests to determine whether the obtained functions represented a signifi cant change in the deviation from a null or a complete model in which the number of parameters was equal to the total number of observations (McCullagh & Nelder, 1989;Dobson, 1999).
A non-parametric Kruskal-Wallis test was used to determine whether there were signifi cant differences both in the total abundance of butterfl ies among the habitats and in the abundance of each butterfl y species among the habitats. For this analysis, habitat type and total abundance of butterfl ies (encompassing both, interactions and simple sightings of individuals) were considered independent and dependent variables, respectively. All analyses were carried out using STATISTICA 10.0 program (StatSoft, 2011).

Recording effort and weather
The analysis of both recording effort curves for butterfl ies and plants indicate that the percentage of species re-corded reached at least 70%, so the data is signifi cantly accurate.
In addition, the statistical relationship between the number of butterfl y interactions and climatic data indicates it is not signifi cantly affected by either daily humidity or temperature, thus there are no climatic biases in the recording of data.

Structure of biological communities
The cluster analysis of vegetation revealed two main fl oristic groups (Fig. 2). Group A included the plots corresponding to the seral and open vegetation. Group B the rest of the plots, mostly corresponding to Q. pyrenaica forest. Within group A, some aggregation patterns were identifi ed. Most of the plots with Nardus meadows, Rubus shrubland and Cytisus shrubland were close to each other.
The butterfl y species richness (both in interactions and sightings) was 43, with 32 species identifi ed per year and 21 (about 40%) species present both years (Table S4). The habitat with the lowest mean number of species in both years was Quercus forest (mean ± standard deviation: 0.50 ± 0.38; signifi cant differences were observed in 2017 in the rest of the habitats; Fig. 3), followed by Nardus meadow (1.44 ± 1.10), Poa grassland (2.56 ± 0.10) and Cytisus shrubland (2.70 ± 0.26). Rubus shrubland had the highest mean butterfl y species richness, but it differed signifi cantly between the two years (8.83 ± 1.17 in 2017; 4.50 ± 2.95 in 2018).
The s imilarity analysis of plots based on butterfl y species occurrences did not show any clustering in any year, except for the Rubus shrubland plots, which all appeared within the same clade (Fig. 4). However, the Nestedness Temperature Calculator analysis revealed that the plots were nested according to their butterfl y composition in both 2017 (P = 6.52*10 -20 ) and 2018 (P = 1.05*10 -4 ). This nested structure of plots confi rmed the existence of a single butterfl y community. The most widely shared species among habitats were also those butterfl ies with the highest number of interactions, such as Melanargia lachesis Hübner, Lycaena alciphron (Rottemburg), Issoria lathonia (L.) and Maniola jurtina (L.). The rarest species in the whole area were Hipparchia semele L., Satyrium ilicis Esper and Pyrgus malvoides (Elwes & Edwards).

Butterfl y-habitat relationships
In 2017, Rubus shrubland was the habitat with the highest butterfl y species richness (21), as well as interactions and sightings (567 and 471, respectively), while Quercus forest had the lowest richness scores (4), interactions and sightings (43 and 120, respectively; Table S4). In 2018, both Cytisus and Rubus shrublands were the habitats with the greatest butterfl y species richness (14 and 12, respectively), interactions (296 and 249) and sightings (17 and 23). The Quercus forest was again the habitat with the lowest scores in butterfl y species richness (3), interactions (8) and sightings (9). Inter-habitat Kruskal-Wallis tests re-vealed signifi cant differences in butterfl y habitat selection in both years (P < 0.0001). A more in-depth species-byspecies Kruskal-Wallis analysis indicated that most butterfl ies showed signifi cant habitat selection in 2017, but not in 2018 (Tables 1 and 2). (Rivas-Martínez et al., 2011). The open (seral) plots with the same type of vegetation were mostly closely grouped by the cluster analysis, although some plots had intermediate fl oristic characteristics, which is common for a natural system with a mosaic of interspersed habitats.

Cluster analysis of vegetation revealed a sharp differentiation between forest (closed habitat) and seral vegetation (open habitat), which is in accordance with the vegetation series described for the Guadarrama Mountains piedmont
In addition, based on our results, distributions of species of butterfl ies were nested at a local landscape scale. The same butterfl y community exploited all habitats, but not equally, throughout the Q. pyrenaica-related vegetation series. Although nested distributional patterns of butterfl ies are usually documented at a regional scale (Fleishman & Murphy, 1999;Franzen & Ranius, 2004; Öckinger et al., Table 1. Summary of the signifi cant relationships between butterfl ies and plants in 2017. For each relationship, the percentage of explained variation is shown and the number of asterisks indicates the P-value: < 0,05, < 0,01 and < 0,001 (! -marginal P-value < 0,1). The signs (+ or -) indicate a positive or negative relationship. ns -non-signifi cant.  Slancarova et al., 2015). A lthough data were recorded at the same period of the year, there was only an overlap of 40% in the species of butterfl y recorded in the two years of this study. Just as the coverage of plants (and, therefore, of nectariferous fl owers) is not the same from one year to the next, even in the same period of the year, the presence of one or another species of the butterfl y community is not the same. There are species whose presence is more stable, since they are more generalists, have a greater dispersal capacity or are more widely distributed throughout all habitats and throughout the year. However, there are other species that are only present when habitat conditions are particularly suitable, as they are more sensitive to various environmental factors, such as weather (Montagud & García Alama, 2010). In 2017, D. deltoides, T. pulegioides and Ca. carpetanus were only visited by two or three species of butterfl y compared to the nine species attracted to R. ulmifolius. However, in 2018, the butterfl y species attracted by these fl owering plants was more equitable (four species for T. pulegioides, Ca. carpetanus and R. ulmifolius, only one for D. deltoides). Our results, however, underline the low number of species of fl owering plants exploited by butterfl ies in the area studied, which is in line with the fi ndings of other similar studies (Jennersten, 1984;Loertscher et al., 1995;Severns et al., 2006;Ezzeddine & Matter, 2008).
The specifi city of butterfl ies with respect to different nectariferous plants has been differently explained, by variations in nectar availability or optimization of feeding (Corbet, 2000); for example, the genus Rubus in another study is reported as a source of an abundance of nectar compared to other species (Holl, 1995). Our study revealed that when R. ulmifolius is scarce, butterfl ies used other plants as an alternative. This indicates that there is a trophic preference gradient, along which butterfl ies optimize foraging by looking for the most optimal cost-benefi t plant. A trophic preference of butterfl ies is reported for a small number of fl owering species, which depends on their availability (Goulson et al., 1997;Lebeau et al., 2017;Szigeti et al., 2018).
A s nectar resources are dispersed throughout an area and concentrated in certain habitats, butterfl ies are not distributed homogeneously in all habitats. This study revealed that open seral habitats, specifi cally woody shrublands, were the areas most widely used by butterfl ies. In those habitats where there was a concentration of resources, butterfl ies also temporarily concentrated. Thus, even though the density of vegetation may be an obstacle to counting butterfl ies it increases the ease of encountering them and their visibility (Nwokwu & Sanderson, 2009).
Rubus shrubland only occupies a tiny area in the area studied, yet its use was disproportionate (Fig. 1). This preference is due to the optimal foraging for R. ulmifolius nectar in a small habitat. When the fl owering of R. ulmifolius ceases, butterfl ies compensate by visiting Cytisus shrublands and Poa grasslands, which contain the other plants they visited. Thus, the distribution of nectariferous plants determines that of butterfl ies at a local scale (Loertscher et al., 1995;Pöyry et al., 2009;van Halder et al., 2011;Fartmann et al., 2013). In some studies, no relationship is reported between plants in fl ower and the distribution of butterfl ies in habitats (Sharp et al., 1974;Holl, 1995). Thus, it is diffi cult to generalize about the diversity of butterfl ies in mosaic habitats based on a short study. However, we believe that butterfl ies gradually change their use of habitats based on the availability of certain fl owering plants.
Cytisus scoparius, which is most characteristic of Q. pyrenaica forest-edges (Rivas-Martínez et al., 2011;Sánchez-Mata et al., 2017), was the most abundant plant in fl ower in 2018, although no signifi cant interaction with butterfl ies was observed. Flowers of Cy. scoparius are not visited by butterfl ies, as they are neither attractive nor accessible since the fl owers do not produce nectar and some force is required to open them (Parker, 1997;Parker et al., 2002).
Mature successional stages (forests and shrublands) are often considered to be poor habitats for butterfl ies (Grill et al., 2005;van Swaay et al., 2006;Kitahara et al., 2007;Colom et al., 2021) and based on the present results the lowest butterfl y richness was recorded in oak forest. Nevertheless, this habitat includes specialized species, such as species of Hipparchia and Pararge (Schmitt, 2003;Bobo et al., 2006). In contrast, Rubus shrubland includes a fl owering plant w ith the highest number of trophic relationships with butterfl ies in early summer, which makes these areas the most interesting habitats within the Guadarramean Mountains in terms of conservation. In addition, as the butterfl y species community was nested across habitats, it can be inferred that protecting Rubus shrubland would be a priority for safeguarding the highest diversity of butterfl y species in the area, followed by Cytisus shrubland and to a lesser extent other grasslands. Thus, apart from oak forest, the woody habitats of seral vegetation has a strong positive effect on butterfl y species richness (Dennis, 2004;Marini et al., 2009).
This study, however, was carried out in a specifi c season and although it was the month when most species and individuals were known to be present, it is necessary to carry out studies in more locations and other seasons to determine the importance of habitats in landscapes. Furthermore, our results indicate that butterfl ies change habitat depending on the resources available in each habitat, and when one habitat is poor in nectar they feed in another. This indicates that habitat patches throughout an area act as complementary resources (Dennis et al., 2003;Oliver et al., 2010) and are a useful buffer against changes in the phenology of nectar resources. Habitats in European mountains are changing due to the loss of traditional practices of fi re, logging and grazing livestock, which results in the homogenization of the landscape in mature forests and the disappearance or reduction of the seral stages of grasslands and thickets (Slancarova et al., 2016;Ubach et al., 2020). This work highlights the importance of heterogeneity and the advantage of a succession of different stages for butterfl ies (Slancarova et al., 2014).

CONCLUSION
For the fi rst time this study compared the butterfl y-fl ower interactions in a mosaic of related habitats in an ecological succession. The same butterfl y community was recorded in all the habitats, but the phenology of the butterfl y-plant interactions changed as the different species of butterfl ies selected habitats based on the availability of nectar from certain plants. Table S1. Average and standard deviation of the temperature (T, in °C) and humidity (H, in %) recorded by the AEMET (State Meteorological Agency, Spain) stations in Colmenar Viejo, Rascafría, P. de Navacerrada on the days in July 2017 and 2018 when the records were collected. The global average and standard deviation (SD) are also shown.  Table S2. Average cover (in %) of the fl owering plants grouped by habitat in 2017. The average total was corrected based on the area of the habitat in the area studied using QGIS.
Year 2017  Table S3. Average cover (in %) of the fl owering plants grouped by habitat in 2018. The average total was corrected based on the area of the habitat in the area studied using QGIS.