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

Study Extent

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 (this dataset). Four major periods in butterfly studies can be recognized: i) an initial study in 1918; ii) short field surveys in the mid–19th 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.

Sampling

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.

Quality Control

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.

Method steps

1. 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.

1. Lepidoptera
   rank: order
2. Hesperiidae
   rank: family
3. Lycaenidae
   rank: family
4. Nymphalidae
   rank: family
5. Papilionidae
   rank: family
6. Pieridae
   rank: family
7. Riodinidae
   rank: family

Galičica is a calcareous mountain range at the junction of Macedonia, Albania and Greece and is shared between the two first mentioned states. 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).