Introduction

Changes of biodiversity patterns along environmental gradients have been the focus of great attention in the context of ongoing global climate changes during the past decades. Short local climatic gradients provide a wide range of habitats suitable for species with different ecological demands and strategies, which increase α-biodiversity in small parts of a landscape. Their marked effect on the structure of soil invertebrate communities has already been documented in numerous groups of meso-/macrofauna, including mites, springtails, myriapods and woodlice1,2,3,4. However, studies focused on latitudinal and elevation gradients often do not provide a sufficient explanation for the distribution and ecological preferences of soil fauna due to their not being concordant results [e.g.,5,6]. Studies dealing with gradients in soil microclimate parameters such as temperature and moisture, mostly associated with pH, C/N ratio, soil texture or nutrient availability appeared more important with adequate explanations for variations in soil faunal biodiversity on a fine spatial scale [e.g.,2,7,8].

Geomorphologically diverse karst landforms such as ravines, valleys, dolines or cave entrances are highly suitable for studies on microclimatic gradients in association with the distribution of the local soil fauna. Cave entrances are suitable objects for observing changes in the local climate because they represent dynamic environments with characteristic microclimatic variations9 and are also complex transition habitats between subterranean and surface environments, supporting high local diversity and endemism. Moreover, these sites have the potential to maintain permanently cold and wet conditions throughout the year, thus serving as refugia for uncommon, specialized species [e.g.,3,10,11,12,13].

The Silická ľadnica Ice Cave is located in southeast part of Slovakia, approximately 1.5 km southwest of Silica village, 12 km south of Rožňava town and 60 km southwest of the eastern Slovak metropolis of Košice (Fig. 1). It is the lowest perennial ice cave located in the northern temperate zone characterized by a noticeable microclimatic inversion in the collapse doline in front of the cave14, which is reflected in the specific distribution of vegetation (Fig. 2). Thanks to the specific microclimate, suitable conditions are provided for psychrophilous plants and animals, which are typical for higher mountain ranges. From the study locality more detailed ecological data on soil fauna are restricted to springtails3,7, while data on pseudoscorpions15, mollusks16, some species of mites17,18 have just faunistic relevance. Therefore, we focused on exact documenting the diversity of oribatid mites as the dominant group of soil mesofauna within this exceptional locality with a very specific character of the local climate, which rarely occurs in a karst landscape of temperate zone.

Fig. 1
figure 1

Location of Silická ľadnica Ice Cave within territory of Slovakia and Slovak Karst National Park. Map sources: Map of Slovakia—https://sk.m.wikipedia.org/wiki/S%C3%BAbor:Slovakia_outline_map.svg; map of Slovak Karst modified according to: https://www.npslovenskykras.sk/wp-content/uploads/2018/01/navstevny_poriadok_mapy.pdf.

Full size image
Fig. 2
figure 2

Study sites location at the Silická ľadnica Ice Cave, Slovakia. Diagram is showing location of the research sites (red circles) in the collapse doline of the ice cave (modified according to Droppa (1962). The gradient went from a scree slope with short vegetation (S1, S2, S3) through a terraced area with shrubs (S4) towards a deciduous forest (S5, S6, S7). For more detailed explanation of study sites, see subsection Site description in Methods).

Full size image

In the present study, we examined for the first time an extensive dataset of Oribatida communities as dominant soil microarthropod group from the Silická ľadnica Ice Cave, originating over a 3-year period. We hypothesized that (1) communities of the model mesofauna group will significantly differ along the microclimatic gradient from the cave entrance slope covered of vegetation towards the adjacent deciduous forest, and (2) the communities at colder sites of the gradient will comprise not only cold-tolerant species from adjacent forest sites but also specialized forms adapted to the harsh conditions at these sites (longer freezing, lower amount of nutrients). The aims of this study were: (1) to reveal the distribution pattern of oribatid mites at seven sites along the microclimatic gradient at the collapse doline of Silická ľadnica Ice Cave; (2) to determine the driving environmental factors shaping the distribution of oribatid communities at study sites (S1–S7); and (3) to assess the function of the studied collapse doline in terms of possible refugium for cold-adapted and endemic Oribatida species.

Results

Soil-chemical and microclimatic data

The soil temperature gradient was recorded with increasing mean, minimal and maximal values from the entrance of the Silická ľadnica Ice Cave up the slope along the transect line (Table 1). Mean soil temperatures showed significant differences between the study sites. Figure 3 shows the soil temperature regime (September 2006–October 2007) at individual sites with similar temperature values and their minimal fluctuations during the late autumn and winter months and, on the other hand, differing values and their pronounced fluctuations (especially at sites S6 and S7) during the warm summer months. Soil moisture increased from the beginning of the transect (with a maximum at site S3) and decreased again towards S7. The chemical characteristics of the soil showed different trends (Table 1). The soil pH/H2O values were mostly neutral at the sites; only at S6 and S7, the soil acidity was lower. The cold sites S1 and S2 had higher C/N ratios than the other sites along the gradient, followed by the wooded sites S5 and S6. Lower values of this parameter were recorded at S3 and S4, with herbal cover, and the lowest in the thermophilous woody site S7.

Table 1 Topographic, vegetation, soil-chemical and microclimatic characteristics of study sites in the Silická ľadnica collapse doline.
Full size table
Fig. 3
figure 3

Mean daily soil temperatures at sites in the collapse doline (9.9.2006–8.9.2007). Explanations: S1–S7: number of study sites.

Full size image

Community abundance

In total, 9610 adult oribatid mites belonging to 142 species and 71 genera from 45 families were collected at all sites within the study period (see also Supplement 2). Juveniles were counted (3962 ind.) without species identification and were thus not included in community analyses with the exception for that of the mean abundance. The Kruskal-Wallis and subsequent Dunn’s post hoc tests confirmed the significant differences (H = 149.4, p < 0.0001) in the abundance of juveniles between particular sites (Fig. 4A). Juvenile mean abundances at the three coldest sites were significantly lower than those from the other sites. Moreover, significantly higher abundances were observed at two warm sites compared to sites in the middle of the gradient. Very similar significant differences (Kruskal-Wallis and subsequent Dunn’s post hoc test, H = 165.2, p < 0.0001) were also recorded in mean abundance of adult oribatids along the study transect (Fig. 4B). Significant differences were also confirmed between adult oribatid mites and juveniles at every site (Mann-Whitney test, p < 0.0001), with abundance means of juveniles being considerably lower. The cold sites (S1–S3) showed not only the lowest oribatid abundances but also a high spatial clustering of individuals, presented by high portion of the samples without oribatid mites in comparison with warmer sites (S4–S7).

Fig. 4
figure 4

Mean abundance per sample at sites in the collapse doline over three years (2005–2007). (A) Juveniles. (B) Adults. Explanations: Different lowercase letters above boxplots indicate sites that significantly differ (p < 0.05) from the others according to Dunn’s post hoc test. S1–S7: number of study sites.

Full size image

Species richness and diversity

Species richness at particular sites ranged from 11 to 74 (Table 2), with a high total number of species (142). The highest richness (74) was detected at one of the warm sites of the gradient (S6). In contrast, the lowest species richness (11) was detected at the cold site (S3). Comparing the mean species richness (species richness/sample) between the sites, the Kruskal-Wallis test showed significant differences (H = 173, p < 0.0001). Dunn’s post hoc test revealed significant differences similar to those recorded in oribatid abundance between cold sites (S1–S3), transition sites in the middle part of gradient (S4, S5) and warm forest sites (S6, S7) (Fig. 5). Mean species richness at S1 site was insignificantly higher than at the other cold sites (S2, S3) with the total number of species at S1 (35) even higher than that of S4 (33) situated in the middle part of the microclimate gradient.

Table 2 Community parameters of adult oribatid mites at sites in the Silická ľadnica collapse doline over three years (2005–2007).
Full size table
Fig. 5
figure 5

Mean number of species (number per sample) at sites in the collapse doline over three years (2005–2007). Explanations: Different lowercase letters above boxplots indicate sites that significantly differ (p < 0.05) from the others according to Dunn’s post hoc test. S1–S7: number of study sites.

Full size image

The diversity indices (H; J’) showed similar trends: higher diversity (H > 2.5) and evenness (J’ ~ 0.90) were recorded at warm, forest sites S6 and S7. In contrast, low diversity (H < 0.3) was noted at the cold, unforested sites S2 and S3 with the most unbalanced community at S2 (J’ = 0.13) (Table 2).

Species distribution along gradient

Dominant species (D > 5%) were Pilogalumna tenuiclava (Berlese, 1908), Conchogneta dalecarlica (Fosslund, 1947), Platynothrus peltifer (C. L. Koch, 1839), Achipteria nitens (Nicolet, 1855), A. coleoptrata (Linnaeus, 1758), Chamobates voigtsi (Oudemans, 1902), representing 46% of total abundance of oribatid adults (see Supplement 3 for illustrative photos of dominant species). The species composition at the cold sites (S1–S3) was highly variable between particular samplings, often with no identical species in two consecutive samplings. Indicator species analysis (IndVal) was used to assess the indicator species, typical for the study sites. For cold sites S1–S3, with very unstable community composition, no indicator species could be specified (Fig. 6). In contrast, Oribatella dudichi Willmann, 1938 was a typical community member at S4, while P. tenuiclava and Ceratoppia bipilis (Hermann, 1804) showed considerably lower relationship with that community. Euzetes globulus (Nicolet, 1855) at S5 and Hermanniella dolosa Grandjean, 1931 at S6 showed the highest percentage relationship to the site community. S7 was the only site with several unequivocal indicator species defined: Multioppia glabra (Mihelcic, 1955), Metabelba propexa (Kulczynski, 1902), Oppiella (Rhinoppia) epilata Miko, 2006, Oppiella (Oppiella) nova (Oudemans, 1902) and Scheloribates (Hemileius) initialis (Berlese, 1908) presented 100% relationship to the site community.

Fig. 6
figure 6

Indicator species analysis (IndVal) based on species presence at sites within every sampling occasion. Explanations: S1–S7 represent study sites. Species abbreviations: ACCA (Achipteria coleoptrata), ACNI (Achipteria nitens), BRBC (Berniniella bicarinata), BRHA (Berniniella hauseri), CDDA (Conchogneta dalecarlica), CHBO (Chamobates borealis), CHVO (Chamobates voigtsi), CRBP (Ceratoppia bipilis), DDRI (Damaeus (Damaeus) riparius), ERHP (Eremaeus hepaticus), EZGO (Euzetes globulus), HEDO (Hermanniella dolosa), HLII (Scheloribates (Hemileius) initialis), HRGB (Hermannia gibba), LIXY (Liacarus xylariae), MEPR (Metabelba propexa), MEPV (Metabelba pulverosa), MOGA (Multioppia glabra), NRAU (Neoribates aurantiacus), OLEP (Oppiella (Rhinoppia) epilata), OLNO (Oppiella (Oppiella) nova), OLOB (Oppiella (Rhinoppia) obsoleta), OLSU (Oppiella (Rhinoppia) subpectinata), ORCA (Oribatella calcarata), ORDU (Oribatella dudichi), PDTE (Pilogalumna tenuiclava), PLPE (Platynothrus peltifer), QOPC (Quadroppia pseudocircumita), SCLV (Scheloribates laevigatus), SCQU (Scheloribates quintus), SGCA (Steganacarus (Tropacarus) carinatus), TCMB (Tectocepheus minor), XETE (Xenillus tegeocranus). The color scale represents the indicator value of species from 0 (blue color) to 100% (red color). The data were Bonferroni corrected to reduce the instance of a false positives. Species with p < 0.05 were boxed as indicator species. The size of the ovals is directly proportional to the percentage values of the species. The larger the oval, the greater the indicative role of the species.

Full size image

Influence of environmental factors on species distribution along study gradient

We examined the relationships between elevation, chemical factors, mean soil temperature and moisture, mean abundance and total number of species of adult mites. Soil moisture negatively correlated with elevation, mean soil temperature (for both parameters: Spearman’s r = − 0.86), and mean abundance of adults and species number of Oribatida (for both last parameters: Spearman’s r = − 0.93). Mean temperature correlated positively with elevation (Spearman’s r = 1). Mean abundance correlated positively with elevation and mean temperature (for both parameters: Spearman’s r = 0.86).

The principal component analysis (PCA) (Fig. 7) performed on data of 33 species in 42 samples (7 sites x 6 sampling dates) showed a high association of species with individual sites. On the other hand, communities were more or less stable in time, with little impact of the year or season. Moreover, the ordination plot depicted the effect of the microclimate gradient: high similarity between communities of (1) cold sites S1 and S2, (2) transition sites S4 and S5, and (3) communities of warm sites S6 and S7. PCA analysis based on pooled data of 33 species for seven sites (7 sites x 1 pooled sampling) showed a similar ordination pattern (Fig. 8A), with separated communities of cold sites (S1–S3) and more similar communities of site couples S4–S5 and S6–S7.

Fig. 7
figure 7

PCA ordination plot of species communities of adult oribatid mites at seven sites (S1-S7) sampled over three years (05–2005, 06–2006, and 07–2007) and two seasons (spring and autumn) in the collapse doline. Variation proportion explained by the first two PC axes (PC1 and PC2) are reported in percentage (in square bracket).

Full size image
Fig. 8
figure 8

Ordination diagram of communities of adult oribatid mites at seven study sites (S1–S7). (A) PCA ordination plot (variation proportion explained by first two PC axes shown in brackets). Vectors of environmental characteristics were passively fitted to ordination plot after PCA. (B) RDA ordination plot with vectors of environmental characteristics explaining constrained variation after adjustment (variation proportion explained by the first two RDA axes and statistical significance shown in brackets). Environmental variables: Tmin—minimum soil temperature, M—soil moisture, pH/H2O and C/N soil ratio. Only subset of species is named to increase clarity. Species abbreviations: ACCA (Achipteria coleoptrata), ACNI (Achipteria nitens), CDDA (Conchogneta dalecarlica), CHBO (Chamobates borealis), CHVO (Chamobates voigtsi), CRBP (Ceratoppia bipilis), DDRI (Damaeus (Damaeus) riparius), EZGO (Euzetes globulus), HRGB (Hermannia gibba), MEPV (Metabelba pulverosa), MOGA (Multioppia glabra), NRAU (Neoribates aurantiacus), ORCA (Oribatella calcarata), ORDU (Oribatella dudichi), PDTE (Pilogalumna tenuiclava), PLPE (Platynothrus peltifer), SCLV (Scheloribates laevigatus).

Full size image

RDA performed on data from 33 species (7 sites × 1 pooled sampling) and environmental data (7 sites × 4 environmental characters) revealed a model, which was statistically significant (p < 0.01). The proportion of explained variation in RDA is relatively high, R2adj = 0.65, which means that most of the total variation in species community composition can be explained by variation in environmental characteristics, and thus RDA was modelled with high confidence. The first two statistically significant RDA axes explained 72.13% of the constrained variation. Three environmental characteristics (pH, M, C/N) were statistically significant (p < 0.01), while Tmin was significant only marginally (p = 0.048). The RDA model showed the same grouping pattern of sites as that in the PCA ordination (Fig. 8B). We thus found support from PCA and RDA ordination pattern clarifying that similarities between the sites in terms of species community composition were considerably driven by environmental factors.

Discussion

Mites are among the most abundant and diverse microarthropods in litter and soil, along with Collembola [e.g.,19]. The distribution pattern of Oribatida communities along the microclimate gradient at the entrance of the Silická ľadnica Ice Cave differed considerably from the pattern observed in Collembola by3 with the highest abundances of springtails and the lowest species richness and diversity indices at the two cold sites of the gradient (S1–S2). It is generally known that habitats with low soil temperatures contain abundant populations of springtails20, but only a few species are adapted to these conditions21. In contrast to springtails, there were no oribatid species bound to this type of microhabitat at sites S1–S3, indicator species analysis did not confirm any indicator species in these cold sites, unlike other sites.

In our study, the abundances and species richness of oribatid mites continually grew along the transect from cold to warm sites. Generally, oribatid communities display lower abundance and species richness in localities with cold microclimate with poor vegetation [e.g.,4,22,23]. Oribatida showed an obvious positive relationship of latitude with species richness, increasing from high latitudes (Antarctic and Arctic regions) towards warm sites of low latitudes of subtropic regions24. However, species richness is positively corelated also with the size of the study area. Nevertheless, environmental temperature is one of the important factors influencing the oribatid species richness. Oribatid mites are considered K-strategists25, since they produce a low number of offspring and have a lower reproduction rate; however, they are stronger competitors in more predictable environments26. Therefore, they are not able to constitute abundant populations at the cold bottom of the Silická ľadnica Ice Cave entrance near the soil permafrost.

Since Oribatida are numerous especially in highly organic soils25, we could assume that the low abundance of the mites at the cold sites (S1–S3) may be associated with the low litter layer and shallow layer of soil. C/N ratio is also considered as factor controlling the organic matter breakdown27. A high C/N ratio over 20–25 at S1 and S2 sites indicates low decomposition rate and low quality of organic substrate in the soil, which could result in low Oribatida abundance at these sites. But C/N ratio values at other sites did not correspond with oribatid abundances, thus this parameter did not explain satisfactorily Oribatid abundance distribution along study gradient.

Various authors have documented that oribatid community composition does not appear to be seasonal28,29. Our three-year study showed that communities did not show significant seasonal changes between spring and autumn in different microclimatic conditions along the gradient, which was shown by PCA ordination. Some studies have pointed out that seasonality can be demonstrated in specific microhabitats. For example30, found non-prominent seasonal changes in abundances of Oribatida in litter, dead wood and mosses, which were explained by the microhabitat instability during seasons. On the other hand31, did not find such recurring patterns in similar microhabitats.

In our study we found that microclimatic factors clearly ruled the distribution pattern of Oribatida across the entrance slope of the cave. In general, soil microclimate has a significant impact on oribatid mite communities19,32,33. We recorded a steep soil temperature gradient in the collapse valley of the Silická ľadnica Ice Cave, from the coldest site located on a steep scree slope close to the floor ice to the warmest site in hornbeam-beech forest on the karstic plateau, with the reverse pattern in the soil moisture. Deep karst valleys and dolines are characteristic with a toposequence similar to the latitudinal zonation of biomes from tundra to deciduous forests [e.g.,3,34], thus representing obvious habitat partitioning across microclimatic gradient at a small space scale.

The coldest site of the gradient (S1) was situated near the perennial ice, within conditions similar to the soils in the alpine zone of mountains in high elevations with cold microclimate and poor vegetation cover, where oribatid mites have low abundances and species richness compared to temperate forests or foothills of elevational gradients4,5,23. They can be also limited in activity or even decimated by low soil temperatures35,36. However, we found a noticeably rich community at this site in comparison with other cold sites (S2, S3), with different values of abundances and species richness across the seasons as well as every separate year. The explanation for this can be found in the location of this site just under the ceiling edge of the cave entrance, thus it was randomly enriched by falling organic materials (leaves, tree branches) from surface of the above karst plateau. For example, the arboreal species Poroliodes farinosus (C. L. Koch, 1840) appeared only at this site, and several species (Damaeus (D.) crispatus (Kulczynski, 1902), D. (D.) riparius Nicolet, 1855, Achipteria nitens, Conchogneta delacarlica, Scheloribates (Hemileius) initialis occurred here in low abundance but clearly preferred warm forest sites (S6, S7). The low abundances of juvenile and adult oribatid mites at the other cold sites (S2, S3) of the gradient confirmed their high sensitivity to low temperatures. At the same time, we found a significantly lower abundance of juveniles compared to adults at all the monitored sites, because mite nymphs are more or less sensitive to lower moisture and temperature levels37,38. Also, juvenile oribatids have longer quiescent pre-molt stage during all juvenile instars which can cause loss of juveniles during extraction39. Less sclerotized juveniles are also attacked more often by predators, e.g., Mesostigmata mites19. Different pattern with decrease of adult density but not the juveniles was also observed40. This could be one of the reasons of low juvenile abundance at sites S2 and S3 where Gamasina abundance was several times higher than Oribatida in all sampling dates. Nevertheless, we consider influence of Gamasina predators on Oribatida density as not the main factor influencing the Oribatida distribution along study transect. There was not obvious relationship between the Gamasina and Oribatida abundances on other study sites. The similarity of the communities at the cold sites was supported by high spatial clustering of individuals and the absence of indicator species. Oribatid mites were considerably more abundant from the middle site S4 of the gradient, with higher species richness due to better microclimatic conditions, well-developed rendzina soil and dense herbal vegetation. It is well known that vegetation determines the organic matter in soil and increases available nutrient content for soil communities41, thus increasing their abundance and diversity36. At this mildly cold and wet site, the montane species Oribatella dudichi was recorded in markedly high abundance. In general, oribatid mites are prone to drought and prefer environments with high soil moisture content33. The hornbeam woody S5 site slightly differed from S4 in ecological conditions and community composition with Euzetes globulus as an indicator species that is known to prefer fresh forest soils42,43. Ordination analyses revealed a higher level of similarity between warm sites S6 and S7 with common forest species, but significantly differing in abundance. Although the soil temperatures at both sites were similar, they differed in the soil parameters that affected oribatid community composition. Hermanniella dolosa was an indicator species at S6 site; it prefers the wetter soils of deciduous forests but is highly resistant to soil drought42,44. The community at this site was made up of species preferring warm, relatively dry habitats with the slightly acid pH of a deciduous forest. In this type of habitat, leaves and relative dead leaves provide a food resource and reduce water loss45. As a result, the community structure was relatively stable despite the high fluctuations in the soil temperature observed at this site.

Oribatid mites represent very rich group of soil microarthropods with a large number of documented species and high species diversity. In this research, we noted the presence of 142 species in the entrance doline of the Silická ľadnica Ice Cave. Some of the specimens were not identified unequivocally to species level because of their damage during manipulation, lack of suitable identification keys or very difficult taxonomic evaluation (genus Phthiracarus). We assumed a priori that no of recorded species is potentially new for the science. Distinct morphological characters of forms unidentified on the species level allowed us to consider them as separate species when calculated species richness and diversity indices. The most recently published list of oribatid mites in Slovakia by46 includes 544 species from the whole country supplemented with several new data18,47,48,49. However, this number will be definitely larger because only 138 species collected within the Silická ľadnica Ice Cave transect are covered by the mentioned lists. Therefore, the remaining four detected species are potentially new to the territory of Slovakia and this must be verified after a thorough study of their morphology. Also some of unclearly identified species, especially Suctobelbella sp.1, S. sp.2, S. sp.3 or Scheloribates aff. laevigatus and Damaeus aff. gracilipes, are probably new for Slovakia. Our findings revealed that this small area hosts relatively high Oribatida diversity which is 26% of the total number of species recorded in Slovakia so far. It may be evidence that karst landforms with pronounced environmental gradients cover a high α-diversity and may represent key habitats for rare species. Specifically, we noticed a Carpathian endemic, Oribatella dudichi, which was also observed in the High Tatra Mountains48 and the Vihorlat Mountains50, mostly in forest soil. Furthermore, the unique status of the Silická ľadnica Ice Cave environment was supported by the occurrence of the troglophilous Pantelozetes cavatica (Kunst, 1962), which is inhabitant of several Western-Carpathian caves. It occupied one of the two coldest sites of the gradient, but surprisingly also one of the warmest sites (S6). This oribatid species is apparently a glacial relict51. Overall, cave entrances represent ecotonal transition habitats between the surface and subterranean environment and form optimal conditions for species preferring one or the other environment [e.g.,9].

More than half of the Oribatida species richness, including endemic and relict species, occurred in the collapse doline of the Silická ľadnica in a low number of individuals (fewer than 10). The constantly changing living conditions due to climate change may cause such small local populations to become extinct. Although soil mites are known to be highly resistant to climate change, especially to soil warming, we can expect certain shifts in soil communities if we consider a long-term period (> 20 years; )38. Šupinský et al.52 studied ice dynamics in the Silická ľadnica Ice Cave and found a significant decrement of the small glacier situated in a deeper section of the cave entrance zone in the recent years. Caves are evidently linked with the surrounding environment, so the loss of ice is a solid indicator of long-term climatic and hydrological changes in the karst landscape. Local habitat heterogeneity, induced by a steep environmental gradient in a karst doline, may play a vital role in providing an indispensable refuge for the soil biota under future climate warming or extreme climatic events38,53. The results provided are intended to be the first step in a long-term study of changes in soil communities on the locality, caused by increasing air temperature and melting ice of small glacier behind the cave entrance. Ordination analyses in our results suggest that soil temperature and moisture have important impact on soil oribatid communities. Changes in these factors caused by e.g., human activities can lead to significant changes in soil fauna communities and, consequently, soil quality. Therefore, cave entrance zones in general should be given attention in regard to conservation priorities. Additionally, many caves are made accessible to tourists, which requires raising awareness and ensuring protection of such unique small-scale environments.

Methods

Site description

The Silická ľadnica Ice Cave (48° 32’ 58.9″ N and 20° 30′ 13.9″ E) is located on the Silická planina Plateau, in the Slovak Karst, which is a part of the Slovak-Aggtelek Karst situated on both sides of Slovakia and Hungary border. This cave with a permanent glacier in its entrance zone is unique phenomenon, being an ice cave situated at the lowest elevation in the temperate zone. The location of the area represents the junction point of the continental and oceanic climate zones. The mean annual air temperature ranges from 5.7 to 8.5 °C and the mean annual atmospheric precipitation from 600 to 700 mm54. This remarkable cave was formed in the Middle Triassic Wetterstein limestones of the Silica Nappe. Its entrance zone is a steeply inclined spacious cavity with a large opening to the surface (corrosive-collapsed abyss). The entrance lies at an elevation of 470 m, facing north towards a thermophilous deciduous forest. The entrance zone, at a distance of approx. 50 m inside the cave, is characteristic by ice decoration developing during the winter and lasting usually until the end of spring. Behind this zone permanent ground ice is deposited in the form of a small glacier, the age of which is at most 2000 years52. The upper edge of the collapse valley-like depression extends to an elevation of 503 m a.s.l14. The entrance zone of the cave is characterized by a strong local climatic gradient analogous to the latitudinal climatic gradient from tundra to deciduous forests, representing clear habitat partitioning on a small spatial scale with diverse habitats for the soil fauna.

At the Silická ľadnica Ice Cave, seven sites (described below) were selected along a 129 m transect line from the bottom of the cave portal to the upper part of the collapse doline (Fig. 2). They were situated on a southwest-oriented scree slope of 5–35° inclination at an elevation of 462.1–499 m a.s.l. with a distinct temperature gradient and a microclimatic inversion. Vegetation at these sites was characterized by Šuvada and Petrášová (unpubl.):

Site 1 (S1) was on a scree slope made up of very stony debris, located 7 m from the edge of soil permafrost—a specific habitat near the cave entrance (see also Supplement 1). This site had pioneer vegetation ass. Cardamino-Chrysosplenietum alternifolii with moss and liverworts and a maximum soil profile depth 3 cm. Site 2 (S2) was on a scree slope with pioneer vegetation, ass. Cardamino-Chrysosplenietum, sparse mosses, liverworts and dense growth of Chrysosplenium alternifolium, and shallow soil 4 cm deep. Site 3 (S3) had pioneer vegetation, ass. Cardamino-Chrysosplenietum alternifolii, sparse growths of mosses and liverworts, herbal cover with Urtica sp. and a maximum soil profile depth of 10 cm. Site 4 (S4) was situated in a microdepression on a moderate slope with well-developed soil 15 cm deep. The site had dense herbal cover, ass. Cardamino-Chrysosplenietum alternifolii and Lunario-Aceretum. Site 5 (S5) was hornbeam-maple wood near a rock wall on a very moderate slope, with dense herbal cover, ass. Lunario-Aceretum and stony soil with a profile to a maximum depth of 3 cm. Site 6 (S6) was a young hornbeam wood with sparse Acer pseudoplatanus trees on a steep slope, sparse herbal cover ass. Waldsteinio-Carpinetum and a soil profile 5–10 cm thick. Site 7 (S7) was a cornel-hornbeam wood on the edge of a karst plateau with moderately dense herbal vegetation ass. Waldsteinio-Carpinetum and a soil profile 15–20 cm thick. During the winter and early spring, sites S1–S3 were partly covered by perennial snow and the soil was frozen.

Field sampling and species determination

Samples were taken once in spring and once in autumn in the years 2005–2007 (18 May and 13 November 2005, 18 May and 25 October 2006, 1 May and 6 October 2007). A total of five soil samples were taken from each site. The samples represented soil cores 10 cm in diameter to a maximum depth of 8 cm. The litter layer was included in soil samples if it was at the site. At the sites with a thin soil layer (S1 and S2), samples were taken with a shovel up to the volume of the extraction cylinder (diameter 10 cm, height 8 cm), so that the volume of all samples was approximately the same. The samples were transferred individually into plastic bags and transported to the laboratory. Soil microarthropods were extracted from the samples in a modified high-gradient Tullgren-type apparatus55 for 7 days and preserved in 75% denaturated ethanol. Acari were sorted under a Leica EZ4 binocular stereomicroscope and identified under a Leica DM1000 light microscope (Leica Microsystems GmbH, Wetzlar, Germany). Oribatid mites were studied on temporary slides using 40% lactic acid to clarify the specimens and determined using common identification keys42,56,57,58. The specimens are deposited in the collection of the Department of Zoology, Institute of Biology and Ecology, Faculty of Science, Pavol Jozef Šafárik University in Košice, Slovakia. The raw dataset generated and analyzed during the current study is available as Supplement 4.

The topographic, soil-chemical and microclimatic characteristics

Table 1 provides data on the topographic, soil-chemical and microclimatic characteristics of the studied sites. Soil temperature was measured every four hours (6 measurements per day) from September 2006 to October 2007 by iButton DS1921G data-loggers exposed at 7–10 cm depth at the given site. The significance of differences between sites was tested using One-way ANOVA after confirmation of normal data distribution using the Shapiro-Wilk normality test. Soil moisture content was analyzed gravimetrically and calculated as an average of three one-time measurements on different dates (October 2007, May 2008 and October 2008). Soil pH was measured potentiometrically using a glass electrode and a reference calomel electrode as active pH in water (pH/H2O). Total soil carbon (C) and nitrogen (N) content were measured by the “wet” method according to59. These three soil-chemical characteristics were measured once in June 2009.

Community data analyses

The following were computed as the main community parameters: A—mean abundance; S—species richness; Smean—mean species richness; H’—Shannon diversity index; J’—Pielou index of evenness and D—dominance (relative abundance). Relationships between microclimatic characteristics and community parameters were calculated using the Spearman’s correlation coefficient and the statistical significance was tested using the t-test. To determinate the significant differences in abundance values of oribatid mites, and the soil temperature and soil moisture between the sites, the Kruskal-Wallis test was applied with Dunn’s multiple comparison post-hoc test. The Mann-Whitney test was used to compare differences between the abundance of adults and juveniles at each site individually. Indicator species analysis (IndVal) was used for identifying species indicative of given groups at the sites. The PAST 4.13 software package60 was used for the aforementioned analyses.

Principal component analysis (PCA) and redundancy analysis (RDA) were used to characterize the associations of oribatid mite communities and effects of environmental factors at the individual sites. Species present in fewer than 50 individuals altogether were excluded from the ordination analyses due to their low predictive value, with 33 species remaining in final datasets. PCA analysis was used to inspect species communities at the sites (S1–S7), between years (2005–2007), and between seasons (spring and autumn sampling). To create a data matrix, species counts were pooled from five samples per site into one, resulting in a final set of 42 samples (7 sites × 3 years × 2 seasons). The final raw species composition data matrix (42 samples × 33 species) was subjected to PCA. PCA was run on a centered Euclidean distance matrix computed on a Hellinger-transformed species abundance matrix in the vegan package (ver. 2.6-461), in R statistical environment (ver. 4.2.262), using the rda() function.

The following environmental data were collected for the sites: elevation (Ele), mean soil moisture (M), minimal soil temperature (Tmin), mean soil temperature (Tmean), maximal soil temperature (Tmax), soil acidity pH (pH) and the C/N ratio (C/N). These data were subjected to correlation analysis in order to exclude those characters, which are highly mutually correlated. A linear Pearson correlation coefficient (r) was used in the calculations. All the environmental characteristics fitted the normality of data (p > 0.05) except Tmin, the data distribution of which slightly deviated from normality (p = 0.038) based on the Shapiro-Wilk test. In subsequent analyses, we kept only those environmental characteristics which were not correlated to each other at |r|=0.75, e.g., M, Tmin, pH and C/N (Ele, Tmean, and Tmean were highly negatively correlated to M). Correlation analyses were performed in the PAST 4.13 software package.

Furthermore, we constructed a data matrix of 7 sites × 33 species (all samples per site were pooled into one) and a corresponding environmental standardised data matrix of 7 sites × 4 environmental characteristics (H, Tmin, pH and C/N). Firstly, we performed PCA to test whether the environmental data may explain differences in species composition data between sites. We used the envfit() function of vegan to overlay the environmental characters passively on the PCA ordination and to test the statistical significance of each of them. Finally, RDA was performed in vegan using the rda() function with default settings to formally test for the ecological interpretation of the species communities on sites. The proportion of explained variation was calculated by using adjusted R2 values (R2adj). Analysis of variance (ANOVA) was applied, and statistical significance was assessed using the permutation test with 999 random permutations.

Conclusions

The microclimatic gradient in the karst doline of the Silická ľadnica Ice Cave is a unique phenomenon providing a wide range of microhabitats within a short distance. We conclude that the local soil microclimate gradient and associated habitat heterogeneity had a considerable effect on the abundance, species richness and distribution pattern of soil Oribatida. The abundances of oribatids were significantly different between the cold and warm part of the gradient. The same pattern was observed in the species richness of adults. PCA and RDA ordination analyses documented a clear delimitation of communities in relation to the soil microclimate and soil-chemical parameters, with most of the recorded species preferring the warm end of the gradient. The studied karst doline revealed unusually high α-alpha diversity of Oribatida with cold and wet sites at its bottom representing key microhabitats (refugia) for rare relict and endemic species. Overall, our study emphasized the importance of the local microclimate for structuring oribatid mite communities and illustrates the need to examine the function of cave entrance zones in maintaining the relict species as unique part of local biodiversity, especially in connection with the long-term climate changes. We are aware that this was a unique ecosystem, and the value of the study itself is not in its scalability, but exactly in documenting and analyzing very rare, or even unique habitats. The presented data may serve as the base for a long-term study of the soil fauna community dynamics at this unique locality, potentially evoked by changes in global climate. These are also the first detailed community data of Oribatida diversity across strong microclimatic gradient at the cave entrance in Central Europe.