of and aspect on diversity on Effects of elevation gradient and aspect on butterfly diversity on Galičica Mountain in the Republic of Macedonia (south–eastern Europe)

Effects of elevation gradient and aspect on butterfly diversity on Galičica Mountain in the Republic of Macedo nia (south – eastern Europe). The patterns of butterfly diversity and community changes in relation to elevation are an interesting and well–covered topic in ecology, but the effects of aspect have rarely been evaluated. Here we studied the changes in butterfly species richness and communities along the elevation gradient and aspect of Galičica Mountain. As expected, species richness changed with altitude, showing a bimodal pattern with two peaks and a declining trend towards higher altitude. Changes were well–correlated with the area in each altitudinal zone, while the effects of productivity were less clear. Butterfly communities at higher altitudes were the most distinct when grouped according to β diversity estimates, followed by mid– and low–altitude communities. Indicator species were found in mid–altitudes and for the combination of low–mid and mid–high altitudes, but not among aspects. Overall, aspect produced a less conclusive effect on species richness and community composition. South and north accounted for most of these differences despite dominant western and eastern and exposition of the mountain slopes. The community temperature index declined with altitude and on the northern aspect, showing these areas hosted more cold–adapted species. Notes on butterfly conservation are provided as 23 species known from historical surveys have not been recorded recently. Resumen


Introduction
It is generally agreed that species diversity declines with altitude, somewhat repeating the more stable latitudinal pattern (Rahbek, 1995). However, changes with elevation are more complex and species diversity tends to follow four main patterns with rising altitude: decreasing, low plateau, low plateau with mid-elevation peak, and mid-elevation peak (McCain et al., 2010). These patterns could be explained by several ecological factors, including climate and productivity, species-area relationship, mid-domain effect, effects of ecotone, biotic factors, evolution, and historical circumstances (Colwell and Lees, 2000;Lomolino, 2001;McCain, 2007;McCain et al., 2010). While studying the effects of biotic and evolutionary factors requires detailed planning and study design, other factors can often be tested more easily. Relatively large mountain ranges with diverse geological structure and climate (from Continental to Mediterranean) make south-eastern Europe an interesting region for studies of elevation gradients in species diversity. At the same time, data on butterfly diversity in this region is scattered in a multitude of faunistic papers, while few studies address the butterfly diversity pattern (Mihoci et al., 2011;Zografou et al., 2014Zografou et al., , 2017Kaltsas et al., 2018). Similarly, the effect of aspect is only rarely taken into consideration when addressing the distribution patterns of butterflies in mountains (Gutiérrez, 1997;Mihoci et al., 2011). However, the aspect of the mountain slope is considered an important topographic factor (i.e. Bennie et al., 2008) and it could play a role in shaping butterfly communities.
The climate along the elevation gradient is known to affect insect biology (Hodkinson, 2005) and is predicted to produce the diversity pattern with midaltitude peaks in temperate regions (Despland et al., 2012;McCain, 2007). Butterflies have shown the same patterns of diversity changes in the mountains as other animal groups, and these changes are commonly guided by variation in climate (MacNally et al., 2003;Gallou et al., 2017). Harsh conditions at a higher elevation have prompted numerous ecological adaptations in butterflies (Junker et al., 2010;Kevan and Shorthouse, 1970;Leingärtner et al., 2014) and it has been suggested that species communities at high elevations tend to show some evolutionary constraints (Pellissier et al., 2013). As the current global climate change induces rapid shifts in butterfly communities across continents (Devictor et al., 2012, Zografou et al., 2014, mountain systems have become increasingly important as refugia for species retracting pole-wards (Fleishman et al., 2000;Hampe and Petit, 2005;Wilson et al., 2007). These rapid changes can be traced by climate change indicators such as the community temperature index (CTI) (Schweiger et al., 2014;Zografou et al., 2014). This index also allows the comparison of extant communities within altitudinal gradient or aspects, where CTI values are expected to decline with altitude and on northern slopes.
We compiled a check-list of butterfly species (Papilionoidea) and a large dataset of all records for the Galičica mountain range in the south-western part of the Republic of Macedonia, combining data from the literature and the authors' field observations. The main goal of the study was to examine how butterfly species richness changes according to altitude and aspect of the mountain, which, in contrast to the nearby ranges, stretches in a north-south direction. We tested the effect of the area, ecosystem productivity and mid-domain effect on elevational patterns of species richness. In addition, we studied the similarities between butterfly communities and determined whether communities along the elevation gradient and aspect differed in community temperature indexes. We also provide notes on the potential extinction of several butterfly species in the area and discuss the conservation value of this mountain range.

Study area
Galičica is a calcareous mountain range at the junction of Macedonia, Albania and Greece and is shared between the two first mentioned states ( fig. 1). It rises between the lakes of Prespa in the east and Ohrid in the west. The climate is mild-continental and Mediterranean, with climatic conditions stabilized by the presence of large water bodies. The mountain is situated on an elevated plain, with the base at 695 m and reaching the altitude of 2,265 m (Avramoski et al., 2010;Ćušterevska, 2016). It stretches in a south-north direction; the east and west are therefore dominant aspects. In 1958 the Macedonian part of Galičica was proclaimed a National Park, with a total area of 24,151 ha, and protected by the state (Matevski et al., 2011).

Data collation and preparation
A historical overview of butterfly fauna of Mt. Galičica was summarized by Krpač et al. (2011). This publication was geo-referenced within Biologer.org biodiversity database (Popović et al., 2020) and used as a baseline for our study. Besides records from the literature, we used an original dataset from Verovnik et al. (2010), unpublished data by I. Jugovic, A. Keymeulen, N. Micevski and the authors' personal records. Species observations of MP, RV and geo-referenced records from the literature can be downloaded after registration on Biologer.org or accessed through GBIF (Doi: 10.15470/jacl7y). Four major periods in butterfly studies can be recognized: i) an initial study in 1918; ii) short field surveys in the mid-19 th century; iii) detailed inventory work in the seventies and eighties; and iv) intensive field studies by several experts from 1995 onwards. A detailed checklist of recorded species is given in supplementary material, differentiating historical (before 1995) and recent data.
Compiling species observations from different datasets could produce bias in the data and should be taken with caution. Thus, before proceeding to the statistical analysis, species occurrence records were subset to those with coordinate precision no less than 2 km. The unique combination of the species name, locality name and year was used to create samples -lists of observed species-and was assumed to be equal to the list of species recorded at a single sampling event. The incidence of frequency of species was determined by combining the samples from a certain elevation zone, or aspect. To remove accidental individual observations, data were further subset to include only samples with more than five observed species.
Altitudes and aspects were extracted from a digital elevation model (European Environment Agency, 2016), while all GIS calculations were made using raster package in R (R Core Team, 2019). Altitude at each point of observation was extracted and used to divide samples in n classes of equal length between maximal and minimal altitude. Calculations were made starting from n = 20 and decreasing the value until good estimates were obtained (i.e. low standard errors, high sample coverage, and good representation of all classes). This resulted in 10 altitudinal classes from 689 m to 2,234 m. Aspect values were calculated in degrees and transformed to four major aspects (cardinal directions): north (0º-45º) and 315º-360º), east (45º-135º), south (135º-225º) and west (225º-315º). To determine dominant mountain aspects for each sample, we included the area in a radius of 1 km around the observation point (i.e. not only the aspect at the exact sampling location). The joint influence of several aspects was minimized by selecting samples with single aspects contributing more than 50 %, and discarding data for samples with more uniform aspect contributions.

Species richness
To compare differences in species diversity between altitudinal classes and aspects, we used species richness measure as a representation of α diversity. Calculations were made using the iNEXT package in R, allowing construction of both rarefaction and extrapolation curves; this provided a robust estimate of the true species diversity even in cases with low and unequal sampling effort (Chao et al., 2014;Hsieh et al., 2016). Incidence frequencies prepared in the previous step were used as input data and are available in supplementary material.
The effect of an area on species richness was assessed by plotting the available area in each altitudinal class or aspect versus estimated species richness. The relationship between productivity and species richness was examined by plotting values of normalized difference vegetation index (NDVI) versus estimated species richness. NDVI was obtained from MODIS (Didan, 2015) between April and September 2017-2019 (representing the vegetation period for the last three years). Average NDVI values were then calculated for each altitudinal class and aspect. Since productivity is directly related to climate, productivity-richness relationship could reflect the influence of climate on butterflies along the altitudinal gradient (Levanoni et al., 2011;Pettorelli et al., 2005). The presence of mid-domain effect was checked by plotting 95 % CI of the null model (1,000 replicates) against observed species richness, using rangemodel R package in R (McCain, 2003;Colwell, 2008). Where applicable, statistical significance was checked using Spearman correlation test in R.

Changes in butterfly communities
To estimate the similarity between butterfly communities (β diversity) at different altitudes and aspects of Mt. Galičica, we used the probability version of the Chao-Jaccard index, provided by CommEcol package in R. This index is less biased than classic similarity indices and it is not sensitive to the omission of some species from the samples (Chao et al., 2005). The results are shown as unrooted dendrograms produced in ape package in R.
In addition to providing the Chao-Jaccard estimate, we used the indicator value index to search for indicator species of a certain elevation zone and aspect (Cáceres and Legendre, 2009). Elevation zones were grouped according to the estimates suggested by Chao-Jaccard index in three elevation classes (low, mid and high altitudes). Estimates were obtained using package indicspecies in R using IndVal.g estimator with 1,000 permutations.
To test whether climate had any effects on shaping the butterfly communities we calculated the community temperature index (CTI) for each delineated altitudinal class and aspect of Mt. Galičica. CTI values were calculated as an average of the species temperature index (STI), accounting also for butterfly abundance (Schweiger et al., 2014). If the climate affected the distribution of butterflies in communities, it could be expected that CTI would decline from the southern to the northern aspect and from lower altitudes towards the top of the mountain. Differences were statistically tested using ANOVA and pairwise t-test in R. Since the test could be affected by unequal sampling, different season, or year of observation, we used only our recent field observations that were collated with an even sampling effort and over the same period.

Results
A total of 4,137 occurrence records from Mt.Galičica were collated, most of these (2,376) being new field observations (see also the dataset published through GBIF (Doi: 10.15470/jacl7y). After removing duplicates and imprecise records, we retained 2,883 observations for the analysis. A total of 168 species were recorded for Mt. Galičica, 159 of which were known before 1995, while 145 were confirmed in recent studies (supplementary material). The decline in species was also evident from an estimated 152 ± 15 species before 1995 to an estimated 143 ± 3 species in the recent period (sample coverage 0.90 and 0.99 respectively).     fig. 2A). We observed a strong correlation between the size of the area available in each altitudinal zone and estimated species richness ( fig. 2A); this relation was statistically significant (rho = 0.81, S = 32, P = 0.008). A comparison of species diversity among different aspects of the mountain showed a more uniform pattern, with somewhat lower estimates for species richness on the northern and eastern slopes and higher richness estimates for the southern and the western slopes ( fig. 2B). The correlation between estimated species richness and the area available in different aspects of the mountain was not significant (rho = -0.20, S = 12, P = 0.917).

Republic of Macedonia
Most of the area of the mountain comprises western and eastern slopes due to the mountain predominant south-north direction ( fig. 1, 2B). However, it should be noted that observed and estimated species richness showed a large discrepancy for the southern aspect, while large standard errors and small sample coverage indicate imprecise calculations. We were unable to determine a clear correlation between ecosystem productivity and butterfly species richness along the altitude, although estimates were close to significant (rho = 0.60, S = 66, P = 0.073) and some positive relation is clearly visible in the graph ( fig. 2C). The highest discrepancies were near the second peak in species richness ( fig. 2C). NDVI tended to increase to about 1,000 meters, then slowly decline towards the highest altitudes. In contrast, the NDVI index showed a similar distribution along the mountain aspects, with somewhat lower values The relationship between species richness at different altitudes and aspects of Mt. Galičica and several predictor variables. Upper graphs show the comparison of available area (histogram) and species richness (lines) for each altitudinal class (A) and aspect (B). The middle graphs depict a similar comparison for productivity (histogram of mean NDVI values) and species richness (lines) between altitudinal class (C) and aspect (D). The lower graph (E) shows the 95 % confidence interval for the null model prediction (blue area) constructed to test the mid-domain effect. Estimated values of species richness are shown as diamonds with lower and upper confidence intervals. Observed species richness values are indicated by circles. Altitudinal classes: a, b,c,152;d,1,307;e,1,462;f,1,616;g,1,770;h,1,925;i,1,080;j,2,234. Fig. 2. La relación entre la riqueza de especies en diferentes altitudes y orientaciones de la montaña Galičica y varias variables predictivas. En los gráficos superiores se muestra la comparación de la superficie disponible (histograma) y la riqueza de especies (líneas) de cada clase de altitud (A) y orientación (B). En los gráficos centrales se muestra una comparación parecida entre la productividad (histograma de los valores medios del índice normalizado diferencial de la vegetación) y la riqueza de especies (líneas) entre clases de altitud C) y orientación (D). En el gráfico inferior (E) se muestra el intervalo de confianza del 95 % para la predicción del modelo nulo (superficie azul) elaborada para comprobar el efecto del dominio medio. Los valores estimados de la riqueza de especies se señalan con rombos con los límites inferior y superior del intervalo de confianza. Los valores de la riqueza de especies observada se indican con círculos. Altitudinal classes: a, 689-844; b, 844-998; c, 998-1,152; d, 1,152-1,307; e, 1,307-1,462; f, 1,462-1,616; g, 1,616-1,770; h, 1,770-1,925; i, 1,925-2,080; j, 2,080-2,234.   2D); this could not be correlated with estimated species richness (rho = -0.80, S = 18, P = 0.333). We found no strong evidence for mid-domain effect on species richness as most points fell outside the confidence interval predicted by the null model ( fig. 2E).

Changes in species communities
The Chao-Jaccard dissimilarity index was generally higher for communities at more distant altitudinal levels, with values ranging from 0.15 to 0.7 (more details available in supplementary material). Communities at the highest altitudes were clearly the most distinctive, but those at mid-elevations also tended to form a separate group in the dendrogram ( fig. 3A). Communities on different aspects showed lower differentiation with a dissimilarity index ranging from 0.1 to 0.36 (supplementary material). Eastern and western aspects hosted practically the same communities, while some differences were evident for northern and especially southern aspects ( fig. 3B).
The indicator value index was computed for 127 and 75 species along the altitude and aspect, respectively. Indicator species were found in the mid-altitude zone and for the combination of low-mid and mid-high zones (table 1), while no indicator species were found for the mountain aspects.   3. A, dendrograma de grupos elaborado a partir del índice de similitud de Chao-Jaccard de la diversidad β entre las comunidades de mariposas de diferentes zonas de altitud; B, orientaciones de la montaña; C, diagramas de caja en los que se muestran los valores medios del índice de temperatura comunitaria de los ensamblajes de mariposas a lo largo del gradiente de altitud, D, diagramas de caja en los que se muestran los ensamblajes de mariposas en las distintas orientaciones de la montaña Galičica. En material suplementatario se proporcionan estadísticas detalladas con resultados de la comparación por pares de valores del índice de temperatura comunitaria. Para las abreviaturas de las clases de altitud, véase la fig

CTI CTI
The community temperature index (CTI) values differed significantly between altitudinal zones (F = 4.904, df = 9, P << 0.001) and also between different aspects (F = 3.558, df = 3, P = 0.014). A declining trend was visible in average CTI values towards the higher altitudes, although the zone around the mid-elevations deviated from the linear decline pattern ( fig. 3C). Statistical significance in pairwise comparisons between altitudinal zones was observed only for more distant altitudinal classes, with and the first two altitudinal zones accounting for most of this significance (details in supplementary material).However, the comparison of communities at different mountain aspects showed more uniform distribution of CTI ( fig. 3D), and only the values on the northern aspect were significantly lower (supplementary material).

Butterfly diversity and its potential decline
Despite the relatively small area studied, butterfly fauna of Mt. Galičica in the Republic of Macedonia is extremely rich, with 168 species recorded. This is comparable to other high biodiversity mountain regions in the Balkan Peninsula, such as Stara Planina in Serbia with 167 species (Langourov, 2019;Popović and Đurić, 2014), Stara Planina in Bulgaria with184 species (Kolev, pers. comm.), Shar Planina with 169 species (Jakšić, 1998;Melovski, 2003), Mt. Olympus with 155 species (Pamperis, 2009), Vitosha with 155 (Beshkov, 2014), Rila with 174,and Pirin with 195 (Kolev,pers. comm.). In addition to high butterfly diversity, Galičica is an important area for butterfly  (van Swaay et al., 2010), and 21 threatened species at a national level (Krpač and Dacremont, 2012). This mountain has also been recognized as a prime butterfly area in Europe (Jakšić, 2003).
The absence of new records for as many as 23 species together with a lower estimated number of species in the recent period is noteworthy, however, and might indicate true extinctions of several taxa. The majority of potentially extinct species include habitat specialists, such as Anthocharis damone, Euchloe penia, Tarucus balkanicus, Pseudochazara amalthea and Spialia phlomidis that were known only from the open rocky habitats, once abundant above Ohrid town. This area in particular (authors pers. observ.) and the entire mountain is becoming overgrown due to abandonment of pasturing, with forest cover increasing on Mt. Galičica from 40 % to 58 %, and pastures declining from 50 % to 24 % in the last decades (Despodovska et al., 2012). Similar patterns have been found throughout Mediterranean Europe, with open grassland butterflies declining due to increased forest cover (Slancarova et al., 2016;Ubach et al., 2020). The lack of new records for woodland species such as Limenitis camilla, Neptis rivularis and Erebia ligea is difficult to explain given the increased forest cover. However, these butterflies are linked to more humid habitats, and recent changes in climate -with prolonged droughts and higher aridity-could cause their decline. The same could be true for the extinction of Erebia oeme, one of seven high alpine species (Varga and Varga-Sipos, 2001) recorded for Mt. Galičica.

Species richness pattern
The estimated richness of butterfly species on the elevation gradient of Galičica Mt. shows a bimodal distribution pattern, with two peaks and a declining trend towards the higher elevations. The decline in species richness between the two peaks was evident even when zones between 1,152 and 1,462 m were grouped together for the analyses (results not shown here), showing that the estimates were not caused by unequal or small sample sizes. Altitudinal patterns for butterfly species richness usually show similar, declining trend towards high elevations with peaks in mid-altitudes (Gutierrez and Menendez, 1995;Gutiérrez, 1997;Wilson et al., 2007;Levanoni et al., 2011;Stefanescu et al., 2011;Despland et al., 2012;Gallou et al., 2017;and the overview in Kaltsas et al., 2018). However, a bimodal pattern is not a rare phenomenon for butterflies peaking at midaltitudes and it has often been recorded in other studies (Gutierrez and Menendez, 1995;Levanoni et al., 2011;Stefanescu et al., 2011, Gallou et al., 2017. There is also considerable evidence for the monotonic decline in butterfly diversity along the elevation gradient (Mihoci et al., 2011;Leingärtner et al., 2014;Kaltsas et al., 2018), but an insignificant trend (Kaltsas et al., 2018) and increase in species richness (Wettstein and Schmid, 1999;Pyrcz et al., 2009) has also been reported. Note that the increase in species richness with altitude is more likely to be caused by specific habitat composition (Wettstein and Schmid, 1999) or by the study of an exclusively montane taxonomic group (Pyrcz et al., 2009).
Species richness was significantly correlated with the area available in each altitudinal class on Galičica Mt., providing strong evidence for species-area relationship. In a similar manner, richness seemed to follow the ecosystem productivity (NDVI) until the prominent peak at about 1,500 m, which is probably caused by the vast area available in this altitudinal zone (coinciding with a large mountain plateau). Numerous studies have shown that butterfly diversity is correlated with precipitation and temperature (Mac Nally et al., 2003;Acharya and Vijayan, 2015). Productivity is derivative of these variables and is predicted to peak at mid altitudes along with species richness of butterflies, but heterogeneity in productivity was also an important predictor (Levanoni et al., 2011). The mid-domain hypothesis suggests that communities developing at high and low altitudes overlap in the mid-range, creating a peak in species richness, as predicted by theoretical models (Colwell and Lees, 2000;McCain et al., 2010). Although we found weak evidence that species richness follows pure predictions of the mid-domain model ( fig. 2E), it is not impossible that some patterns in species diversity could be caused by overlapping of the communities from higher and lower elevation zones (see the discussion below). Diversity of butterflies along the elevation gradient could also be affected by heterogeneity of habitats, human disturbance (Lien, 2013;Gallou et al., 2017), biotic interactions, and evolutionary history (Pellissier et al., 2013). With so many factors involved, it is more likely that the actual species richness pattern is a product of these factors (Lomolino, 2001).
The effects of aspect on species communities are better studied in plants (i.e. Gallardo-Cruz et al., 2009;Holland and Steyn, 1975) and little is known about its effect on butterfly diversity. On the similarly oriented Toiyabe range in Nevada, eastern slopes are shown to host more butterfly species (Fleishman et al., 1997). This was explained by more diverse habitats on the eastern slopes and better connectivity, allowing intrusion of southern faunal elements. On Mt. Biokovo in Croatia, butterfly richness was higher on northern slopes than on the very steep and vegetation-poor southern aspects, which are also more exposed to strong bora winds (Mihoci et al., 2011). Our results did not provide any conclusive evidence on the effects of aspect on butterfly richness, despite the unequal distribution of the aspects in Mt Galičica. Somewhat lower estimates of species richness were observed in the northern and eastern exposures, but the estimates for the southern aspect were unreliable and discrepancy between observed and estimated species richness was strong ( fig. 2B, 2D).

Community changes
Changes in butterfly communities along the altitude ( fig. 3A) coincided well with the transition zones of major forest communities on Galičica Mt., with oak forest up to 1,200/1,400 m and beech forests up to 1,900 m (Matevski et al., 2011). Interestingly, only the mid-elevation zone (1,462-1,770 m) had unique indicator species (15), and it shared some indicator species with lower (18) and higher elevation zones (5). This mid-elevation zone also had the highest species richness, and at least partially, this richness could be attributed to the overlapping species from higher and lower elevation zones, providing some evidence for mid-domain or ecotone effects. The changes observed in the community temperature index over the altitudes ( fig. 3C) matched our theoretical assumptions. As altitude increased, the communities were composed of more cold adapted species, resulting in lower values of CTI, but these changes were subtle.
In contrast to clear separation of high altitude species by the Chao-Jaccard index, no indicator species were found for butterfly communities in the highest elevation zone as high altitude species were shared with mid-elevations (table 1). Butterfly species at higher elevations are known to have distinctive ecological adaptations that restrict their distribution (Kevan and Shorthouse, 1970;Leingärtner et al., 2014) and some of them are known to have evolved as separate lineages in the Balkan refugia during the last glaciations (Varga and Varga-Sipos, 2001;Schmitt and Varga, 2012). Admittedly, the number of high alpine species is low for Mt. Galičica, but such a pattern is also observed on other calcareous mountains in the Mediterranean region of south-east Europe (i.e. Sijarić, 1983;Mihoci et al., 2011;Kaltsas et al., 2018) and could probably be explained by predominant xeric conditions and limited alpine areas above 2,000 m.
Finally, a small sample size at the highest altitudes in our study and the low number of high altitude species could limit the number of estimable indicator species.
The most striking result regarding aspect of Galičica Mt. is the lack of differentiation in species assemblages between western and eastern slopes even though they are the dominant aspects on the mountain and clearly separated by a central mountain plateau. A similar result was obtained in a study comparing tree cover on Mt. Galičica; it showed that only southern and northern slopes hosted distinctive communities, while western and eastern aspects were similar in species composition (Matevski et al., 2011). This is in line with lower estimates of CTI for butterflies on the northern aspect. Northern slopes receive less insolation (Geiger et al., 1995) and thus host more cold-adapted species. The butterfly assemblage on the northern aspect are therefore likely affected by local ecological factors such as woodland cover, which is more extensive at colder and more humid parts of this calcareous massif.

Conclusion
Altitude was found to be a strong driver in shaping butterfly species diversity, with species richness peaking twice at mid-altitudes and declining towards the top of the mountain. The explanation for this pattern is likely linked to a highly significant species-area relationship, with possible effects of ecosystem productivity (surrogate for climate) and community overlap. According to the β diversity estimates, most distinctive butterfly fauna was found close to the top of the mountain. However, indicator species were confined to the mid-elevation zone or shared between mid-low and mid-high elevations.
The effect of aspect was not as strong and easy to interpret as the effect of altitude. Comparing all the evidence, it can be concluded that eastern and western aspect have similar richness and species composition and some subtle differences were found when contrasting southern and northern slopes, the latter with marginally lower butterfly richness.
Habitat changes on Mt. Galičica in recent decades, especially the decline of pasturing, have already played a strong role in shaping butterfly communities and have probably caused extinction of several habitat specialist butterfly species. Habitats could easily be restored and sustained with traditional grazing, creating a system that would benefit the local communities, maintain the diverse calcareous grasslands, and enable long-term survival of the unique butterfly communities of Mt. Galičica. statistical inference. Ecology, 90 (12)

Supplementary material
Details of recorded species and results of statistical analyses. The dataset contains a checklist of butterfly species recorded for Mt. Galičica, given separately for historical and recent periods, altitudinal classes, and four major aspects of the mountain. The original matrix used to calculate the species richness index is provided and detailed results of statistical analyses mentioned in the paper are given.