Environmental niche divergence of species from Merodon ruficornis group (Diptera: Syrphidae)

In this paper we analyzed environmental differentiation of closely related species from the Merodon ruficornis group. By applying principal component analysis (PCA) and environmental niche modelling (ENM)-based techniques, we estimated the level of niche divergence of closely related species. Our results indicate that ecology has an important role in the diversification process in related species from the M. ruficornis group. Distribution patterns of all analyzed species are mainly affected by the limiting effects of the temperature of the coldest quarter and month, as well as by the precipitation of the wettest and driest quarters. Our results demonstrated that among all related species, with the exception of M. ovaloides, overall or partial divergence in environmental space is present. Importantly, the results indicate that the environmental niches of all endemic species are restricted to smaller parts of the environmental space. In the case of niche overlap, the niches of endemic species are placed along the border of the realized niche of the widespread related species. For species in which distribution is not limited by geographical barriers, environmental preferences could be considered as limiting factors for further expansion, as in the case of M. alexandri, a lowland species with very strict climatic adaptations. Knowledge about the environmental factors that might influence the diversification process can provide an explanation for the high diversity in other Merodon species groups.


INTRODUCTION
Insects are sensitive to environmental changes and respond to all climatic oscillations within their niches.When insect species are confronted with severe environmental changes (e.g.large-scale climatic changes in geological history), they have three possible solutions: evolve, change range or become extinct [1].There are a large number of examples showing insect species that have changed range and/or evolved as a response to repeated isolation during glaciation-interglaciation cycles [2][3][4][5][6][7].In the Western Palearctic, there are multiple glacial refugia: southern Europe (Iberian, Apennine and Balkan peninsulas), the Anatolian Peninsula and southern Caucasus [2,8,9].Fossil remains of different animal and plant species provide a clear picture that during the last ice age organisms with the northern and central European distribution survived in southern refuges beside the Mediterranean Sea [2], on the Iberian, Apennine and Balkan peninsulas, and perhaps some near the Caucasus and Caspian Sea [10].
The genus Merodon Meigen (Diptera: Syrphidae: Merodontini) is comprised of more than 170 species and is the largest hoverfly genus in Europe [11,12].Species from this genus are connected with bulbforming plants which they use as food sources for the larvae [13][14][15], and as a result some are considered as pests [16].This genus is distributed over the Palaearctic and Ethiopian regions [17] and the Balkan and Anatolian peninsulas, which represent the main centers of biodiversity with a large number of local endemic species on islands, mountains or other isolated areas [11,12,18].Species from the Merodon ruficornis group have been the subject of several integrated studies [17,19,[20][21][22][23].In the last decade, several new species from this group were described [18,[20][21][22][23].Recently, this monophyletic group was revised and identification key and distributional data for a total of 18 species were provided.According to the diagnosis given in Vujić et al. [21], species within the Merodon ruficornis group can be divided into several clusters of morphologically closest species.The distribution pattern of each cluster is similar, comprising of one widely distributed species and another sister species with a restricted range.A speciation process has taken place in these geographically isolated populations and has resulted in the formation of closely related species in different parts of the range of the common ancestor, as hypothesized by Vujić et al. [21].
Analyses of niche relationships among closely related taxa can provide insights into the ecological distinctiveness and mechanisms responsible for the diversification [24][25][26][27].
A species' fundamental niche consists of the set of all conditions that allow for its long-term survival, whereas its realized niche is that subset of the fundamental niche that it actually occupies.By definition, then, environmental conditions at the occurrence localities constitute samples from the realized niche.Analysis of species occurrences enables an approximation of the species' realized niche, in the study area and environmental dimensions being considered [28].
Here we used principal component analysis (PCA) and ecological niche modelling (ENM)-based techniques to examine the mechanisms of ecological divergence among closely related species within the Merodon ruficornis group.
PCA of environmental parameters is a method widely used in studying niche determination and differentiation.It allows comparison of environmental niches among investigated taxa to determine if they differ and if they do, to what degree.Additionally, it can isolate environmental parameters correlated with speciation [29] and can contribute to investigation of the environmental limits of taxa, suitable environments for protected or invasive species and species delimitation [7,[30][31][32][33][34][35][36][37][38].
Species distribution models (SDMs) are increasingly being used to address a diverse range of applied and theoretical questions [39][40][41].Also known as ecological niche models and bioclimatic envelopes, SDMs are correlative models that use environmental and/or geographic information to explain observed patterns of species occurrences [42].
The aim of this paper was to i) define the environmental envelope for each investigated species based on distributional, climate and elevation data; ii) quantify and compare environmental niches among related species using three different approaches; iii) establish abiotic factors that influence species range, and iv) identify factors that can correlate with speciation.Results obtained in this paper can contribute to understanding the speciation process in related species and also can explain why the Mediterranean basin is a hotspot for particular hoverfly genera.

Clusters of morphologically closest species
Environmental preferences were identified and compared among morphologically related species that are aggregated into four clusters based on morphological characters according to Vujić et al. [21].Every cluster follows the same pattern of distribution that is characterized by one widely distributed and two endemic species with a narrow range.

Data preparation
Only records with precise distributional data were used.Localities with geographic coordinates were used without modification.Records without geographic coordinates were georeferenced and visually checked using Google Earth (Google Inc, California, USA, https://www.google.com/earth;accessed on Feb 10, 2015) based on detailed locality information.Localities with unclear information were removed from the data set.

Environmental data extraction
Nineteen bioclimatic variables and elevation data for each location were generated on the basis of World-Clim dataset [43] in 2.5 arc-minutes resolution using DivaGis software (DIVA-GIS version 7.5).The bioclimatic variables were derived from the monthly temperature and rainfall values in order to generate more biologically meaningful variables.They represent annual seasonality trends of temperature and precipitation, and extreme or limiting environmental factors.Distribution and richness maps were created in DivaGis.

Point based analysis of environmental variables
PCA was used to define the climatic profiles of the investigated species and to evaluate whether morphologically close species were separating along certain environmental factors.PCA was carried out applying a normal varimax rotation of factor loadings.Only factors with an eigenvalue greater than 1 were considered significant.Climate variables with a factor loading greater than ±0.8 were interpreted as meaningfully correlated with the factor.A scatter plot of PCA score values was used to graphically display the position of the analyzed species in environmental space.ANOVA was carried out to determine overall differences in the derived factor scores between species.Significant difference in factor scores between species pairs was tested using Fisher's LSD post hoc.Indication of the adaptability of species for each PC was represented as standard deviation of factor scores.All statistical analyses were computed in Statistica® for Windows (version 12, Statsoft Inc, Tulsa, OK).

Ecological niche modelling
Niche projections were conducted using MAXENT (version 3.3.3),which is one of the most commonly used algorithms for species distribution modelling [39,[44][45][46][47][48].The maximum entropy (MAXENT) model [44] originates from the statistical mechanics [49].The idea of MAXENT is to estimate the target distribution by finding the distribution of maximum entropy (i.e., that is closest to uniform) [39] by using occurrence records of the species and environmental variables.To eliminate highly correlated variables, we used VIF (variance inflation factor) analysis integrated in the package usdm (R package, version 1.1-15) in R platform (R Core Team, Vienna, Austria), which resulted in different combinations of Bioclim variables for every species used for modelling procedures.70% of data were used for training the model, and 30% were used as test data.Equal weight was given to training sensitivity and specificity.Accuracy of the model was assessed using the AUC (area under the receiver-operator curve) value, which is one of the most commonly used measures of model performance [45].Values of AUC ranged from 0.5 for models with no predictive ability, to 1.0 for models giving perfect predictions [50].

Tests of niche overlap
Niche overlap, identity and similarity tests for each species were calculated using ENMTools [51,52].Niche overlap analyses (ranging from 0 − no overlap to 1 − complete overlap) were applied within each cluster of species using Schoener'D index [53].The niche identity test was used to determine whether the ENMs generated for the species pairs in each cluster were identical or exhibited statistically significant difference (p≤0.05,p≤0.01).We compared the niche overlap value (D) of pairs of Merodon species to a null distribution, with an overlap value=100.The null hypothesis of niche identity was rejected when the empirically value for D was significantly distinct from the value expected from the pseudoreplicated datasets [52].The identity test is conservative as it only assesses if the pair of species tolerates the exact same set of environmental conditions and it does not consider the surrounding space.For this reason, we also used the niche similarity test to evaluate if the examined species were more or less similar than expected by chance, based on the environmental differences in their ranges [51,52].The null hypothesis of niche similarity was rejected if the true measured overlap values were significantly lower (or higher) than the values generated by the background test.This test is conducted in both directions, and different directions may yield different results.

RESULTS
In this paper, diversity studies (species richness and altitudinal zonation) of the M. ruficornis group were analyzed, and for each cluster, environmental niche divergences were measured using three different approaches: i) point based analysis (PCA), ii) niche projections (MAXENT) and iii) niche overlap, identity and similarity tests (ENMTools).
The distribution of species from the Merodon ruficornis group encompasses Western Palearctic, without representatives on the Iberian Peninsula and northern Africa (Fig. 1A).The species altitudinal range was from sea level to approximately 3000 m a.s.l.(Fig. 2A).Widely distributed species, such as M. auripes, M. loewi, M. ruficornis and M. trebevicensis, reside in areas at different altitudes.Endemics and species with a narrow range had smaller altitudinal ranges and were connected with mountain areas, all except M. alexandri and M. gallicus, which are lowland species (Fig. 2A).The relation between altitude and species richness is depicted in Fig. 2B.The highest number of species can be found at altitudes from 1000 to 2000 m.The species richness and diversity map show that the Balkan Peninsula is the main hot spot for this species group, and together with Anatolian Peninsula represent the biodiversity center of the M. ruficornis group (Fig. 1B).

st cluster of morphologically close species
Merodon ruficornis has a wide altitudinal range (0-2000 m a.s.l.), unlike M. abruzzensis and M. lamellatus, which have considerably smaller altitudinal ranges, connected with mountain habitats.According to the available data, M. abruzzensis occurs at 1000-2000 m, while M. lamellatus has a slightly lower altitudinal range, but inhabits higher altitudes, from 1800 to 2500 m a.s.l.(Fig. 2A).
To identify the specific climate variables that best explain the geographic ranges occupied by each species, PCA analysis was used.PCA extracted five environmental dimensions (principal components), which together explained 96% of the variance.Overall significant separation among species occurs along all five PCs (ANOVA: PC1: F 2,201 =30.95, p<0.00000;PC2: F 2,201 =147.39,p<0.00000;PC3: F 2,201 =13.22,p<0.000004;PC4: F 2,201 =9.40; p<0.000125;PC5: F 2,201 =31.20; p<0.00000).The first principal component explains 50% of total variation, and is negatively correlated with altitude data and positively correlated with temperature (Bio1, Bio6, Bio11).Environmental variability connected with temperature variables is explained by PC1, PC2 and PC4, while precipitation was connected to PC3 and PC5.PC1 depicted a gradient for temperature of the coldest quarter and month, but an opposite gradient for the altitude value (Fig. 3A).The second PC was related to the mean diurnal range (Bio2) and isothermality (Bio3), explaining 17% of the total environmental variation.Additionally, PC4 illustrated the temperature gradient in the driest and the wettest quarter with 11% of the total variation (Fig. 3B).PC3 (13%) and PC5 (6%) elucidate precipitation in the wettest quarter and month, and precipitation seasonality, respectively (Table S1).Scatter plots showed species environmental niches and a clear separation of M. lamellatus and partial overlap of M. ruficornis and M. abruzzensis niches (Fig. 3).In order to confirm that there were substantial environmental differences among species pairs, the Fisher LSD test was conducted.Results showed that PC1 and PC2 differentiated all three species (p<0.05).Additionally, M. abruzzensis and M. lamellatus differ from M. ruficornis in PC3 and PC4, while M. lamellatus differs from M. abruzzensis and M. ruficornis in PC5.The climatic adaptability of each species was calculated as standard deviation around the centroids of the PCs (Table S2).Standard deviation for PC1 showed that all three species have relatively similar adaptability to temperature and altitudinal changes.The high value for PC2 revealed a broad adaptability of M. abruzzensis to monthly and daily temperature ranges, unlike M. lamellatus.Merodon ruficornis had a higher adaptability related to temperature in the driest and wettest quarter and precipitation, reflecting its wide distribution (Table S2).
The predicted distribution of M. ruficornis showed that the most suitable environment with a probability of occurrence above 70% is in the Dinaric mountain range, but also in the part of the Alps between Italy, France and Switzerland, in the high parts of the Apennines in Italy and on the mountains in Turkey, Georgia and Russia along the coastal zone of the Black Sea.50-70% of probability was in the Dinaric Mountains at lower altitudes, the lower slopes of the Apennines, central France, central and southern Germany and the lower slopes of the mountains around the Black Sea.Distribution of M. lamellatus is completely separated from the predicted distributions of M. ruficornis, while M. abruzzensis is distributed in the area where the model showed a 50% possibility of occurrence for M. ruficornis, indicating that the ecological niches of these two species are more similar (Fig. 4A).
Based on ENM, niche overlap is higher between M. ruficornis and M. abruzzensis (D=0.442)than between M. abruzzensis and M. lamellatus (D=0.100) and between M. ruficornis and M. lamellatus (D=0.162)(Table S3).According to the identity test, the null hypothesis of niche identity is rejected, mean-  ing that the climate envelopes of all pairwise comparisons of Merodon species are highly significantly more different than expected to randomly occur (TableS3).On the other hand, in the analysis of niche similarity, the null hypothesis is not rejected in all pairs of species.The environmental niches of M. ruficornis and M. abruzzensis were not statistically different (p≥0.05),presenting 44% of geographic overlap.The ecological niche of M. abruzzensis was less similar than expected by chance to that of M. lamellatus in both directions (Table S3).Other pairs of Merodon taxa shared niches that were more similar than expected by chance, but only in one direction (M.lamellatus with M. ruficornis) and not vice versa.

nd cluster of morphologically close species
Merodon auripes is a species with a mostly continental distribution and wide altitudinal range (Fig. 2A).Based on the available material, M. alexandri is a lowland species, while M. ponticus is a high mountain species (Fig. 2A).
The environmental niches of species from the investigated cluster differed in four dimensions (PC1-PC4) and altogether explained 91% of the variability.The results of ANOVA showed that environmental preferences varied significantly across species in three PCs (PC1: F 2,145 =38.77801, p<0.00000;PC2: F 2,145 =222.2665,p<0.00000;PC3: F 2,145 =16.84435, p<0.00000; PC4: F 2,145 =0.246920, p<0.781531).Niche axis 1 (PC1) allowed a clear distinction between all three species.PC1 depicted a gradient of precipitation levels of the driest quarter and month and annual precipitation, and explained 46% of the total variation (Fig. 5A).PC2 and PC3 were connected to temperature variables and explained 24% and 14% of variability, respectively.PC2 separated the lowland species M. alexandri from M. auripes and M. ponticus according to the mean temperature in the coldest quarter and month (Fig. 5A).PC3 clearly separated M. ponticus from M. alexandri and M. auripes based on altitude and mean temperature of the warmest quarter and month (Fig. 5B).Fisher's LSD test showed that PC1 differentiated all three species, while PC2 significantly separated M. alexandri from M. ponticus and M. auripes, and PC3 separated M. ponticus from M. auripes and M. alexandri (p<0.05).These environmental distinctions across environmental axes are visible on a scatter plot (Fig. 5).
Merodon alexandri and M. ponticus have low deviations for each PC, reflecting their narrow geographic distributions and very strict climatic adaptations, especially for M. alexandri (Table S2).Merodon auripes had a wide adaptability to climatic fluctuations in the driest and warmer quarters (PC1, PC3), and the lowest values for temperature fluctuations in the coldest quarters.Based on niche projection, the highest probability of occurrence (70-100%) for M. auripes was along the coastal zone of the Black Sea, on the Carpatho-Balkan mountain belt, Dinaric mountain range, the Apennines and in the mountain area between Czech Republic, Germany and Austria.A probability between of 50-70% was found in the regions that border the areas with the highest percent- age of suitability, and have lower altitudes compared with them.The potential distribution of M. auripes is clearly separated from the distribution of M. alexandri and M. ponticus (Fig. 4B).
In the second cluster, the niche overlap between the investigated Merodon species had generally low levels, from 0.177 for M. auripes-M.alexandri to 0.329 for M. auripes-M.ponticus (Table S3).Randomization tests of niche identity indicated that the species in each pair were more different than expected, so they are not ecologically equivalent (p≤0.01).The niche similarity was less pronounced than expected by chance in both directions for the pair M. alexandri and M. ponticus.For some pairs of Merodon species, niche similarity was greater than expected by random also in both directions, such as M. auripes with M. alexandri (Table S3).No significant differences (p≥0.05) were recorded between M. auripes and M. ponticus in both directions.

rd cluster of morphologically close species
Merodon trebevicensis has a wide altitudinal range from sea level to approximately 2000 m a.s.l.Merodon gallicus has a smaller altitudinal range, from sea level up to 500 m.Merodon hoplitis is a mountain species (300-1450 m a.s.l.).
Taking into account only the factors with eigenvalues superior or equal to 1, four principal components were retained.These four factors together accounted for 90% of total variance.Environmental preferences differ significantly across the species in all four PCs using ANOVA (PC1: F 2,257 =271.56,p<0.00000;PC2: F 2,257 =72.78, p<0.00000;PC3: F 2,257 =9.91, p<0.000071; PC4: F 2,257 =72.10; p<0.00000).PC1, with 52% of total variation, explained the gradient of annual precipitation and precipitation levels in the wettest and driest quarter and month.This axis separated M. hoplitis from M. gallicus and M. trebevicensis.PC2 showed temperature seasonality and annual temperature range.This axis was responsible for 21% of the total environmental variation (Fig. 6A).The temperature gradient of the coldest quarter and month and annual temperature levels were described with PC3 and explained 10% of total variation.This axis separated M. gallicus from the other two species.Precipitation seasonality was responsible for a small part of the total variation (PC4: 7%) and distinguished M. gallicus from M. hoplitis (Fig. 6B) (Table S1).Fisher's LSD test showed that all species pairs differed significantly among the four PCs (p<0.05).Scatter plots showed the position of the investigated species in environmental space, and it was clear that the environmental niches differed partially and complexly, with a small overlap on all four axes (Fig. 6).
Of all three species from this cluster, M. trebevicensis had the widest adaptability.The very high value for PC3 reflects its wide distributional and altitudinal ranges.According to standard deviation values, M. gallicus had high adaptability for precipitation seasonality, and relatively high adaptability for annual and seasonal temperature variation (Table S2).
Potential distribution for Merodon trebevicensis revealed that the highest probability of occurrence (70-100%) was on the Dinaric and Carpatho-Balkan Mountains, in Italy and on the mountains along the Black Sea, and probability from 50-70% was in the surrounding regions at lower altitudes.The potential distribution of M. trebevicensis (70-85%) overlaps with the distribution of M. hoplitis, while distribution of M. gallicus is completely separated (Fig. 4C).
The results of ENMTools showed that the niche overlap of the investigated species in the third cluster were between 0.346 (M.hoplitis-M.trebevicensis) and 0.543 (M.gallicus-M.trebevicensis) (Table S3).According to the identity test, the climate envelopes of all pairwise comparisons of Merodon species were highly significantly more different than expected as a random occurrence (Table S3).Results of the background test also supported ecological differentiation between species pairs; in most of the investigated pairs, niche similarity was greater than expected by chance, except between M. hoplitis and M. gallicus, where there was no significant difference (p≥0.05).

th cluster of morphologically close species
Merodon loewi has the greatest altitudinal range of all species from the group, from 0 to above 2500 m a.s.l., while M. ovaloides and M. turcicus inhabit mountains above 1000 m (Fig. 2A).Point-based analysis of environmental variables showed that environmental preferences between these species differed in four dimensions (PC1-PC4).They together explained 92% of the total variation.ANOVA showed that overall differences across species were significant only in PC1 (F 2,311 =17.92, p<0.00000).Fisher's LSD test showed that environmental preferences significantly differed only between M. turcicus and M. loewi in PC1 and PC2 (p<0.05), and only these axes were included in the results' interpretation.PC1 was responsible for 36% of the total variation and depicted a gradient of annual mean temperature (Bio1), temperature in the coldest month (Bio6) and quarter (Bio11), as well as the mean temperature of the warmest quarter (Bio10) and the altitudinal range (Alt).PC2 is related to the precipitation level in the wettest quarter (Bio16) and month (Bio13), annual precipitation level (Bio12) and precipitation level in the coldest quarter (Bio19) (Fig. 7) (Table S1).PC2 explained 27% of the total variation.The environmental niche of M. turcicus partially differentiated from that of M. loewi, and had a much narrower range of temperature and precipitation than M. loewi.
The high values of standard deviation of PC1 and PC2 revealed the high adaptability of M. loewi, which is in concordance with its distributional and altitudinal range (Table S2).Merodon ovaloides had the lowest adaptability, while M. turcicus showed relatively high adaptation to annual temperature fluctuations (Table S2).
Based on niche projection, the highest probability of occurrence for M. loewi (70-100%) was along the Pindos mountains in Greece, and continued along Mt.Taygetos on the Peloponnese peninsula.The same probability of occurrence was in Turkey, along the Pontic Mountains and Mt.Taurus.50-70% probability was on lower areas around these mountains, also on the Apennines and along the coastal zone of the Adriatic and Black seas.The potential distribution of M. loewi overlaps with the realized distributions of the two other species from this cluster (M.turcicus and M. ovaloides) in the range of 70-85%, which points to the similarities in ecological preferences among this species (Fig. 4D).
Cluster 4 had the highest value of niche overlap between species compared to the other investigated clusters, from 0.431 for M. loewi-M.turcicus to 0.550 for M. ovaloides-M.turcicus.The identity test showed that the pairs of species M. loewi-M.turcicus and M. ovaloides-M.loewi were significantly more different than expected by chance (Table S3), while there were no significant differences (p≥0.05) between M. ovaloides-M.turcicus.For all pairs of Merodon species in this cluster, the niche similarity was more or less similar than expected by chance, but only in one direction (Table S3).

DISCUSSION
Recent research has clearly demonstrated the role of ecology in the speciation process [7,31,37,54,55].It is clear that a species' range is influenced by its ecological niche [56], which is defined by the combination of environmental conditions and resources that are necessary for an organism to maintain a viable population [57].
This paper illustrates a trend of environmental divergence among related species from the Merodon ruficornis group.By observing the positions of environmental space of all investigated species, it was clear that species with a narrow range were restricted to a smaller part of the environmental space and constituted a subset of the realized niche of the widespread species.Low standard deviation values for endemic species reflect their narrow geographic distribution and very strict climatic adaptations, whereas large values for widely distributed species reflect their very broad adaptability.According to this, those widespread species (M.auripes, M. loewi, M. ruficornis and M. trebevicensis) could be considered as possible ancestors.
In all four comparisons our results showed that all species pairs differed in environmental niches, except M. ovaloides from M. turcicus and M. loewi.Merodon turcicus and M. ovaloides have narrow distributions connected to mountain peaks that are part of the potential distribution of M. loewi (Fig. 4D).The niche overlap between M. ovaloides and M. turcicus was the highest compared to the other investigated species, and presented the large degree of geographical overlap, while M. ovaloides and M. loewi do not differ on any PC axes.For all other species pairs, the identity test showed statistically significant niche diversification (*p≤0.05,**p≤0.01),and PCA showed separation on at least one environmental axis.All investigated species differed as regards temperature in the coldest quarter and month, while precipitation of driest quarter and month was responsible for the environmental separation of species from the second and third clusters.Precipitation in the wettest quarter and month separated M. ruficornis from M. abruzzensis and M. lamellatus and all species from the third cluster.These factors represent the extreme or limiting environmental factors and they could be considered as partly responsible for speciation.
The axis related to altitude is always connected with temperature and depicts species replacement along altitude and temperature ranges (Cluster 1: PC1, Cluster 2: PC3, Cluster 4: PC1).Eight of the 12 species can be found at altitudes from 1000 to 2000 m a.s.l.The number of species above 2000 m gradually decreases and mostly encompasses endemic species.Their habitats are usually connected to high mountains in southern Europe, northern and eastern Turkey.Habitats in these areas with varied topography and lower latitude represent likely places for multiple glacial refugia [2,[8][9][10]58,59].Almost every mountain peak in this region has its own endemic species from the Merodon ruficornis group because these areas act like islands, which limits gene flow among them.
During the ice ages, many thermophilic organisms disappeared from most of their previous distribution areas and only survived in areas at lower latitudes with suitable climatic conditions.In Europe, these areas are located south of the transversal European high mountain systems, and represent the three most important refugial areas of the Mediterranean peninsulas [10,58,60,61].Furthermore, the Maghreb and Asia Minor were also amongst the important refuges for thermophilic species during the glacial periods [8,62].It can be assumed that species from the Merodon ruficornis group, like many other investigated taxa [4,10], have undergone many range contractions and expansions in and out of refugia in the south.During interglacials, some populations would survive in the high mountains, and during the cold periods (glaciation) they would descend to populate refugia at lower altitudes [61].
The mountains of the Alps and the Pyrenees acted as significant geographical barriers to the further spread of M. gallicus, while for M. abruzzensis, M. alexandri and M. hoplitis climate may explain the distributional constraints.These three species have significantly different climatic preferences restricted to a small part of the environmental space and low climatic adaptability.Although M. gallicus and M. trebevicensis are extremely morphologically similar, differences in their environmental niches are statistically significant.Merodon gallicus inhabits drier and warmer areas with less seasonal temperature variability, which is not suitable for morphologically related species.These environmental factors, in combination with the geological history of the Maritime Alps (France) could explain the present-day distribution of M. gallicus [63][64][65].Merodon hoplitis is distributed along the coastline of Montenegro and Croatia, as the consequence of newly acquired specific climatic conditions (higher humidity and milder temperatures during the year) that are provided in Dinaric microrefugia near the coastal zone [66].Also, the distribution of M. alexandri is limited by specific environmental factors.This species has very low adaptability to all environmental factors represented by the PCs.It is a lowland species that inhabits Ukraine and the Russian steppes, which are characterized by extreme continental conditions.Species diversity in the Merodon ruficornis group is high in comparison with other groups of this genus.Twelve of 18 (66.5%)species are endemic and have very local distribution connected with only a few mountain peaks on the Balkan Peninsula, Turkey and in the Caucasus region [21].These areas were influenced by diverse climatic and other geographic conditions in the past.This can explain a large number of local endemic insect species in the main refugia of the Western Palearctic, which are also the main hotspots for many insect species [4,60].
Future field investigations are necessary to enlarge the sampled material, particularly for rare species, as well as to complement this diverse group with new members, probably endemics from isolated mountain peaks.As glacial and interglacial periods have left distinctive marks on the genomes of many species, including those from the Merodon ruficornis group, further investigations should utilize recent advances in molecular analysis in order to improve our understanding of the phylogenetic relationships.

Fig. 2 .
Fig. 2. Results of altitude analysis of species from the Merodon ruficornis group.A -variability plot of species altitudinal gradient; B -relation between altitude and species richness.

Fig. 3 .
Fig. 3. Distribution of species from cluster 1 in environmental space.A -temperature and altitude related axes (PC1 and PC2); B -temperature and precipitation related axes (PC3 and PC4).

Fig. 7 .
Fig. 7. Distribution of species from cluster 4 in environmental space.
The first cluster of morphologically close species contained Merodon ruficornis Meigen, 1822, Merodon abruzzensis van der Goot, 1969 and Merodon lamellatus Vujić & Radenković, 2012.Merodon ruficornis is predominantly distributed in central parts of Europe, including France to the west and the Apennine and Balkan peninsulas to the south.Merodon abruzzensis is locally endemic with a range on the Abruzzo Mountains in Italy, at the border of M. ruficornis distribution, while M. lamellatus is an endemic species from the Turkish part of the Caucasian region.Merodon trebevicensis Strobl, 1900, Merodon gallicus Vujić & Radenković, 2012 and Merodon hoplitis Hurkmans, 2012.Merodon trebevicensis is a widely distributed species from the Alps in the west, to Crimea and central Turkey in the east.