Karst landforms as microrefugia for soil Collembola: Open versus forested dolines

Surficial enclosed depressions with a bowl or funnel shape, also called dolines, solution dolines, or sinkholes, are among the most numerous types of surface landforms in karst areas around the world. They usually range from meters to tens of meters in diameter, and their inner slopes vary from subhorizontal to nearly vertical (Sauro, 2019). Such dolines are characterized by the presence of distinctive environmental conditions, which vary considerably from slope to slope and from top to bottom, often resulting in steep microclimatic gradients. The bottoms of the dolines tend to receive more water due to the runoff of precipitation, may retain snow cover longer than the other doline microhabitats and are also a place of cold-air pooling. Furthermore, bottoms and poleward-facing slopes receive less insolation than other parts of the doline and surrounding area, thus providing relatively cool and humid conditions, which are most pronounced at midlatitudes (Holland and Steyn, 1975; Bárány-Kevei, 1999; Kiss et al., 2020).

Recently, numerous studies have documented that such high microclimatic and habitat heterogeneity represents an important factor ruling the local plant (Hortal et al., 2009; Stein et al., 2014; Bátori et al., 2017; Bátori et al., 2021; Yáñez-Espinosa et al., 2022; Bátori et al., 2023a; Bátori et al., 2023b) and animal communities (Báldi, 2008; Seibold et al., 2016; Bátori et al., 2019b; Bátori et al., 2020; Marcin et al., 2021; Bátori et al., 2022; Marcin et al., 2022) in a karst landscape. It has been shown that karst dolines can shape community structure, the distribution of functional/ecological groups of plants and animals (according to temperature and moisture preferences), and provide shelter for specific species that are rare or fully absent from the surrounding plateau (Sólymos et al., 2009; Vilisics et al., 2011; Kemencei et al., 2014; Bátori et al., 2019b; Bátori et al., 2022). A similar effect of microclimatic gradients was documented in springtails (Collembola) as good model organisms for surveying the functional biodiversity of soils (Hopkin, 1997; Rusek, 1998; Potapov et al., 2020) and a reliable bioindicative group for local and regional climatic variations in soil environments (e.g., Lindberg et al., 2002; Kardol et al., 2011). The effect on diversity and distribution of functional/ecological groups of Collembola at microclimatically distinct sites was documented in open/grassland dolines (Marcin et al., 2022), which have more pronounced differences in microclimatic conditions than in forested areas due to the absence of a tree canopy. In contrast, forest cover markedly lowers daily and seasonal temperature variations, resulting in lower soil temperatures and overall soil microclimatic stability compared to open habitats (Wallwork, 1970; Coleman and Crosley, 1995; Bátori et al., 2011; Bátori et al., 2014; De Frenne et al., 2019; De Frenne et al., 2021; Bátori et al., 2022), although it still provides a high diversity of animals and plants (e.g., Sólymos et al., 2009; Vilisics et al., 2011; Kemencei et al., 2014; Maclean et al., 2015; Bátori et al., 2019b; Bátori et al., 2022).

Although the microrefugia character of various types of dolines has been documented (e.g., Raschmanová et al., 2008; Růžička et al., 2016; Raschmanová et al., 2018; Papáč et al., 2019; Marcin et al., 2021; Marcin et al., 2022), there is still a dearth of investigation into the solution dolines and the mechanisms that drive the distributions of animal taxa within these specific habitats in karst landscapes (e.g., Nagy and Sólymos, 2002; Vilisics et al., 2011; Kemencei et al., 2014; Bátori et al., 2019b; Bátori et al., 2022). Various landforms associated with diverse environmental conditions with apparent microclimatic gradients and habitat heterogeneity may be considered natural “habitat islands,” indispensable habitats for various endemic and relict species (Raschmanová et al., 2008; Vilisics et al., 2011; Bátori et al., 2014; Keppel et al., 2018; Raschmanová et al., 2018; Bátori et al., 2019b; Öztürk and Savran, 2020). Such vulnerable species are highly threatened by ongoing climate change. In addition, natural microclimatic gradients in/within dolines may enhance the probability of speciation due to adaptation to diverse environmental conditions. Different microclimates in habitats may create different selective pressures, leading to isolated local populations (Raschmanová et al., 2017). Therefore, dolines that present a dominant feature in karst landscapes have crucial conservation importance. Overall, higher diversity and abundance of Collembola is expected in forested dolines compared to open habitats (Perez et al., 2013) due to the favorable soil microclimate, well-developed soil profile, and presence of forest litter, which provides heterogeneous microhabitats and direct or indirect food resource for soil fauna (Wallwork, 1970; Hopkin, 1997; Gaston, 2000; Rusek, 2001; Bardgett, 2002; Salamon et al., 2004). Moreover, karst dolines are often inhabited by specialized cold-adapted species, most of them considered unique climatic relicts, which may significantly contribute to soil biodiversity. We expect their predominance in forested dolines, where climatic fluctuations are diminished, specifically in cold and wet doline bottoms and N-facing slopes (Kemencei et al., 2014; Bátori et al., 2019b; Marcin et al., 2022).

In this study, we aimed at (1) a comparison of the diversity and structure of soil Collembola communities at sites along microclimatic gradients in karst open and forested dolines, (2) the detection of the responses of individual ecological groups of species regarding their affinities to the soil microclimate, and (3) an assessment of the function of the 2 types of dolines as potential climate refugia for psychrophilic forms of soil Collembola.

The studied dolines are located in the southwestern part of the Silická plateau in the Slovak Karst (Slovakia), a part of the Slovak–Aggtelek Karst geomorphological unit situated in Slovakia and Hungary (Figure 1A). This area has a typical karst character with various landforms, such as gorges, solution and collapse dolines, pits and sinkholes, with a specific microclimatic regime characterized by a more or less developed temperature inversion. The mean annual air temperature in the area ranges from 5.7°C to 8.5°C, and the mean annual precipitation ranges from 630 to 990 mm (Rozložník et al., 1994).

Two open and two forested dolines near the village of Kečovo in the Slovak Karst, situated at an elevation of 415–465 m a.s.l., were selected for the study (Figure 1B–D). All the dolines have a circular shape, while the forested dolines are slightly larger than the open dolines. Open dolines (1) and (2) had a diameter of around 175 and 130 m and a depth of 17 and 7 m, respectively, and forested dolines (3) and (4) had a diameter of around 300 and 250 m and a depth of 35 and 20 m, respectively. The distance between unforested dolines (1) and (2) and forested dolines (3) and (4) was approximately 15 m. The distance between the open and forested dolines was approximately 350 m. The open dolines are situated in an area actively used for cattle grazing, while the forested dolines are overgrown with at least 70-year-old oak forest.

Four sites were selected across each doline from south to north: N—north-facing slope, B—doline bottom, S—south-facing slope, and PB—edge of the doline on the karst plateau. Transects across the dolines differed in length depending on the dimension of the individual doline. The sampling and soil fauna extraction design followed Marcin et al. (2022). A total of 5 soil samples were taken from each site (80 in total) on November 16, 2021. The soil samples represented soil cores 10 cm in diameter, taken from a maximum depth of 8–10 cm (depending on the soil thickness). Collembola individuals were sorted under a binocular Leica S6E stereomicroscope and identified under a Carl Zeiss Axiolab A1 phase-contrast microscope (Carl Zeiss Microscopy, Oberkochen, Germany) to the species level using multiple taxonomic keys (e.g., Gisin, 1960; Stach, 1960, 1963; Zimdars and Dunger, 1994; Jordana et al., 1997; Fjellberg, 1998; Bretfeld, 1999; Potapov, 2001; Thibaud et al., 2004; Jordana, 2012).

The Collembola specimens are deposited in the collection of the Department of Zoology, Institute of Biology and Ecology, Pavol Jozef Šafárik University in Košice, Košice, Slovakia. The research adhered to the conditions of License No. 2661/2017-6.3. from the Ministry of the Environment of the Slovak Republic.

The coordinates of the locations in the dolines within the study sites were recorded in the national coordinate system (S-JTSK; EPSG: 5514) by the global navigation satellite system using a Topcon HYPER HR receiver with a connection to the reference network of the Slovak real-time positioning service in real-time kinematic mode. To provide information on the topography of the individual sites in the dolines, data on the elevation, slope, exposure, topographic index, solar radiation, and insolation were recorded. The topographic index is directly related to the likelihood of water accumulation at a given site (Quinn et al., 1991), and it is used as a proxy for soil moisture estimates in the present study.

In this study, the parameter of solar radiation is defined as the total annual solar radiation received at the soil surface, and insolation represents the annual number of hours of direct solar radiation at a given site in clear-sky conditions. Geomorphometric parameters were derived from the digital terrain model (DTM) using ArcGIS Pro (Spatial Analyst toolbox) and GRASS GIS (r.slope aspect module) software. As the computational inputs, 1-m resolution (DTM) and the digital surface model (DSM—vegetation cover included) were used. Digital models were derived by the linear interpolation of the classified point cloud from airborne laser scanning (provided by the Geodesy, Cartography, and Cadastre Authority of the Slovak Republic) with a declared horizontal accuracy of 0.08 m and a vertical accuracy of 0.09 m. Solar radiation was modeled using the Area Solar Radiation tool in ArcGIS Pro software, with outputs of total radiation and insolation in clear-sky conditions for the whole year for the DSM inputs.

Vegetation associations at the sites were characterized according to (Hegedüšová-Vantarová and Škodová, 2014): Doline (1), the plateau site associated with the Onobrychido viciifoliae-Brometum erecti Müller 1966 association, the N-facing slope associated with the Scabioso ochroleucae-Brachypodietum pinnati Klika 1933 association, the S-facing slope with thermophilous vegetation, the Rosetum gallicae Kaiser 1926 association, and the bottom contained relatively high and dense growths of grasses with characteristic cold vegetation of the Alchemillo-Arrhenatheretum elatioris Sougnez and Limbourg 1963 association; Doline (2), the plateau site associated with the Festuco rupicolae-Caricetum humilis Klika 1939 association, the N-facing slope associated with the Onobrychido viciifoliae-Brometum erecti Müller 1966 association, the S-facing slope with the Scabioso ochroleucae-Brachypodietum pinnati Klika 1933 association, and the bottom with the Pastinaco sativae-Arrhenatheretum elatioris Passarge 1964 association. The forested dolines were characterized by a thermophilous oak wood with hygrophilous herbal vegetation at the doline bottoms (R Šuvada, personal communication, 31/05/2022).

Soil temperature at the sites was measured continually every 4 h by data loggers exposed to 3 cm of soil depth from June 1, 2021, to May 31, 2022. For each site, the annual mean (T_mean), minimum (T_min), and maximum (T_max) soil temperatures were calculated. The differences in soil temperature (T_mean) between the sites for each doline were tested using a Friedman analysis of variance (ANOVA) test.

The soil-chemical parameters were analyzed at each doline site (Soil Science and Conservation Research Institute, Bratislava, Slovakia). Soil pH_H2O was measured potentiometrically using a glass electrode and a reference calomel electrode as an active pH in water (ISO 10390:2005). Organic carbon (COX) and total nitrogen (NTOT) content (ISO 10694:1995) were measured on a CHNS-O Elemental Euro EA 3000 analyser (Italy) according to Kobza et al. (2011).

In order to characterize the Collembola communities at the doline sites, several principal community parameters were calculated: mean abundance and species richness as quantitative parameters, and Shannon diversity and the Pielou equitability index as qualitative parameters. The differences in abundance mean and species richness between individual sites were tested using the Kruskal–Wallis and post hoc test (Statsoft, 2013) for each doline separately. The relationships between the environmental factors and community parameters were estimated for each doline separately using the nonparametric Spearman’s correlation coefficient (Statsoft, 2013). Species with a dominance D >5% were considered numerically dominant (Tischler, 1955).

Communities between the dolines and sites were compared using nonmetric multidimensional scaling (NMDS) ordination based on species abundance at the sites. Autopilot with slow and thorough mode and Sörensen (Bray–Curtis) distance (recommended for community data) were selected. After randomization runs, a 2-dimensional solution was accepted as optimal. The NMDS analysis was performed by the PC-ORD 7 package (McCune et al., 2002; McCune et al., 2016). Species present in less than 3 individuals in the whole material were excluded from this analysis due to their low explanatory value.

Based on the experience of the authors and literature data on the ecology of individual species, Collembola were differentiated into 6 ecological groups according to their habitat preferences in relation to (1) soil moisture, that is, hygrophilous, mesophilous, and xerophilous/xeroresistant, and (2) temperature, that is, eurythermic, cold-adapted, and thermophilous species (Appendix). The community patterns of the ecological groups (species numbers and abundance means) were tested using general linear model analysis, and 5 graphs were used to illustrate the distribution of the ecological groups at the different sites; specifically, 2 models for cold-adapted and thermophilous species were built for temperature, and 3 models for hygrophilous, mesophilous, and xerophilous species for moisture. The model for eurythermic species was not prepared due to its low informative value in terms of the research aims. In these models, the abundance and species numbers of the ecological groups were treated as a dependent variable, the site as a fixed factor, and the location (dolines 1–4) as a random factor. For every ecological group, the significance of the differences between the sites was tested using the Fisher LSD test (Statsoft, 2013).

The values for the topographic, soil temperature, and edaphic parameters at the sites are summarized in Table 1. In the forested dolines, differences in temperature between individual sites were lower than in the open dolines. The Friedman ANOVA test confirmed the significant differences in the soil mean temperatures between the doline sites: (1) χ² = 1,006.6, N = 365, df = 3, p < 0.00001, (2) χ² = 691.8, N = 365, df = 3, p < 0.00001, (3) χ² = 862.1, N = 365, df = 3, p < 0.00001, and (4) χ² = 562.9, N = 365, df = 3, p < 0.00001, respectively.

In each doline, the N-facing slopes and bottoms showed lower soil temperatures than S-facing slopes and plateau sites (Table 1). More specifically, the lowest soil temperature means in open dolines (1) and (2) were recorded at the bottoms and N-facing slopes, respectively, while within the forested dolines these values were lowest at the bottoms. The highest mean temperatures were documented at the S-facing slopes within the open dolines and plateau sites within the forested dolines. The T_max values were much lower in the forested dolines. Markedly less variable soil temperatures were observed at sites in forested compared with the open dolines. The warm sites of the open dolines (1S, 1PB, 2S, and 2PB) in particular showed the highest temperature variations (Figure 2). Warm plateau sites of all the dolines were characterized by the highest values of C_OX and N_TOT, except for doline (1), with the highest values of both parameters on the warm S-facing slope (1S). The overall highest values of these parameters were recorded at the forested plateau sites (3PB and 4PB). Soil acidity (pH) was lower at the sites of the forested dolines compared with the open dolines, showing a decreasing pattern toward the bottoms of the dolines, except for doline (4).

A total of 11,642 individuals from 88 species of Collembola were recorded at the studied sites (Appendix). The abundance means of the communities at the sites across the dolines varied considerably, 10,777–18,803 ind.m⁻² in doline (1), 14,421–21,682 ind.m⁻² in (2), 16,128–47,032 ind.m⁻² in (3), and 7,873–25,096 ind.m⁻² in (4). The number of species in the dolines ranged from 21 to 28 (1), from 18 to 22 (2), from 27 to 39 (3), and from 30 to 34 (4) (Table 2).

Overall, 47 species were recorded in the open dolines. Eighteen species occurred exclusively within the open dolines (Appendix); 17 species occurred exclusively at sites inside the open dolines (N, B, and S), and 2 species occurred only at the plateau sites (PB). In the forested dolines, 69 species were documented in total. Altogether, 40 species occurred exclusively within the forested dolines, 25 species occurred exclusively inside the forested dolines, and 5 species occurred exclusively at the karst plateau.

Eleven Carpathian/Western Carpathian endemics were documented (Appendix), two of them—Jevania weinerae (site 2PB) and Orthonychiurus rectopapillatus (site 1N)—exclusively in the open dolines and other 9 exclusively in the forested dolines, with the largest representation at plateau sites (3PB and 4PB): Deuteraphorura jitkae, Deutonura stachi, Endonura tatricola, Hymenaphorura pseudosibirica, Micranurida bescidica, Micranurida vontoernei, Superodontella tyverica, Tetracanthella montana, and Pumilinura loksai.

Four cold-adapted species were documented in the open dolines, with the overall highest abundance at the bottoms (1B and 2B), and 14 cold-adapted species were found in the forested dolines with the highest abundances at the N-facing slopes (3N and 4N; Appendix). In open doline (1), the community parameters showed an increasing pattern toward the colder sites (1N and 1B), with the abundance and species richness being lowest at the plateau (1PB) and diversity indices at the warm S-facing slope (1S; Table 2). In doline (2), the community parameters showed the opposite pattern, increasing toward the plateau (2PB), except for the abundance mean, which was highest at the bottom (2B). The abundance and the species richness were lowest at the N-facing slope (2N), while diversity indices were lowest at the bottom (2B). Regarding the forested dolines, in doline (3), the community parameters showed an increasing pattern toward the colder sites (3 N—abundance and species richness and 3B—diversity indices), with the lowest values at the plateau (3 PB), except for the equitability index, which was lowest at the N-facing slope (3N). Isotomiella minor, Folsomia manolachei, and Folsomia penicula showed aggregation patterns at site 3N (Appendix), with considerably high mean abundances in a few soil samples with a high proportion of juveniles. Similarly, in doline (4), the community parameters showed an increasing pattern toward the colder sites (4N and 4B), except for the species richness being highest at the plateau (4PB) and lowest at the N-facing slope (4N). An exceptionally low abundance was documented at the doline bottom (4B).

Significant differences in the abundance mean between individual sites were confirmed by Kruskal–Wallis ANOVA only for the forested dolines, in doline (3), H (3, N = 20) = 9.63, p = 0.022, with significant differences confirmed by post hoc test between sites 3N and 3PB (p = .028), and in doline (4), H (3, N = 20) = 8.81, p = 0.032. Significant differences in species richness between individual sites were confirmed only for open doline (2), H (3, N = 20) = 8.10, p = 0.044.

Several significant correlations were revealed between the environmental factors and community parameters. In doline (1), the Shannon diversity index correlated negatively with T_min (R = −1.00, p < .05); the abundance of Parisotoma notabilis correlated negatively with T_min, T_max, C_OX, and N_TOT (R = −1.00, p < 0.05), and the abundance of Lepidocyrtus cyaneus correlated positively with soil pH (R = 1.00, p < 0.05). In doline (2), species richness correlated positively with T_min (R = 1.00, p < 0.05); the mean abundance correlated negatively with soil pH (R = −1.00, p < 0.05); the abundance of Hypogastrura assimilis correlated positively with T_max (R = 1.00, p < 0.05), and the abundance of Protaphorura pannonica correlated positively with C_OX and N_TOT (R = 1.00, p < 0.05). In doline (3), the species richness correlated negatively with T_min (R = −1.00, p < 0.05); the Shannon diversity index correlated negatively with C_OX (R = −1.00, p < 0.05); the abundance of F. manolachei correlated negatively with T_min (R = −1.00, p < 0.05); the abundance of Megalothorax willemi correlated positively with T_min (R = 1.00, p < 0.05), and the abundance of P. notabilis correlated negatively with C_OX (R = −1.00, p < 0.05). No significant correlations between the environmental factors and community parameters were revealed in doline (4).

Collembola communities of the dolines were further analyzed using NMDS ordination. The best 2-dimensional solution had a final stress of 9.73, p < 0.00001 after 54 iterations, which was confirmed by a Monte Carlo permutation test with p = 0.004 and a mean stress of 10.9 for real data and 250 runs for both real and randomized data. The first and second axes explained 73% and 16% of the variance. The NMDS diagram outlined 2 clearly delimited clusters concerning the type of karst doline (Figure 3). Within the open dolines (cluster on the left side of the diagram), a noticeable delimitation concerning the microclimatic character of sites was observed. Warmer sites (S, PB) on the left side of the cluster, with characteristic species, such as Isotoma anglicana, Proisotomodes bipunctatus, and Protaphorura serbica, were separated from the cooler sites (N, B) with characteristic species, such as Lepidocyrtus lignorum, P. notabilis, and Megalothorax minimus. Site 1S was markedly separated from the other open doline sites, presumably due to the exclusive presence of characteristic xerothermophilous species P. serbica and P. bipunctatus at this site. The forested dolines (cluster on the right side of the diagram) did not show any marked delimitation of species concerning doline sites. Cold site 3N, markedly separated from the other forested dolines sites, was associated with abundant F. penicula. Cold site 3B, again markedly separated from the other forest sites, was associated with eurythermic and hygrophilous Folsomia quadrioculata.

In the open dolines, significant differences in the species numbers were confirmed for the cold-adapted, thermophilous, and xerophilous/xeroresistant groups, and regarding mean abundances, significant differences were confirmed for the thermophilous and xerophilous/xeroresistant groups (Table 3). Within the forested dolines, significant differences in the species numbers were confirmed for the cold-adapted, hygrophilous, and xerophilous/xeroresistant groups, and regarding mean abundances, significant differences were confirmed for all ecological groups (Table 3).

Graphs of the number of species occurrences and mean abundance (Figures 4 and 5) showed several characteristic distributional patterns of ecological groups of Collembola at topographically and microclimatically different sites. In the open dolines (Figure 4) regarding species numbers, the cold-adapted group showed a preference for bottoms (B) but mainly for plateau sites (PB), which was significantly different from sites at the slopes of the dolines (N, S). Regarding mean abundances, this group showed a preference for cold doline bottoms (B), although this relationship was not significant. Regarding both species numbers and mean abundances, the thermophilous and xerophilous/xeroresistant groups showed a clear preference for the warm sites, specifically S-facing slopes (S) and plateau sites (PB), which were significantly different from the cold sites (N, B). In the hygrophilous group, no significant differences were observed in species richness and mean abundances; the hygrophilous species did not occur at S-facing slopes. Regarding species richness and mean abundance, the mesophilous group showed the same distributional pattern, with a preference for cold sites, such as doline bottoms and N-facing slopes, but they were nonsignificant.

Within the forested dolines (Figure 5), regarding both species numbers and mean abundances, the cold-adapted group showed a strong preference for cold N-facing slopes, which differed significantly from the other sites. In contrast, xerophilous/xeroresistant species showed a preference for the warmer sites, specifically S-facing slopes (S) and plateau sites (PB), significantly different from the other sites in species richness and also mean abundance. The same pattern was observed in the thermophilous group, which also showed a preference for warmer sites (S, PB) and was significantly different from other sites only in mean abundance. Moreover, hygrophilous species had a significant relationship with cold and wet bottoms of the dolines (B). Similarly, the mesophilous group showed a certain preference for colder and wetter sites (N-facing slopes), although this relationship was not significant.

The results of this study correspond with previous studies showing that relatively shallow surficial karst dolines provide a broad variety of microhabitats (Sólymos et al., 2009; Vilisics et al., 2011; Kemencei et al., 2014; Bátori et al., 2019b; Bárány-Kevei et al., 2021). The relatively small depth of the studied dolines resulted in rather less pronounced microclimatic gradients. Regardless, significant differences in the soil mean temperatures between the sites for all dolines were documented. The differences between the forested sites were considerably lower and also had lower temperature fluctuations compared to the sites of the open dolines.

Moreover, site topography in a doline is associated with the amount of soil nutrients. In our study, plateaus were characterized by high amounts of carbon and nitrogen more prominently in the forested dolines, with a decreasing trend toward the bottoms of the dolines. Soils of steep slopes contain less carbon and nitrogen due to the lower accumulation of organic matter compared to the flat plateaus (Laurance et al., 1999). The low soil carbon and nitrogen content at the bottoms of the dolines may be attributed to nutrient leaching, which increases with humidity (Brouwer and Powell, 1998; Quintero-Ruiz et al., 2019).

Microclimate and habitat heterogeneity play a significant role in determining biodiversity patterns in terrestrial fauna in karst dolines (Báldi, 2008; Sólymos et al., 2009; Vilisics et al., 2011; Kemencei et al., 2014; Seibold et al., 2016; Bátori et al., 2019b; Marcin et al., 2021; Bátori et al., 2022). However, the diversity and abundance patterns across karst dolines cannot be uniformly applied to all terrestrial invertebrates due to their group-specific and especially species-specific relations to the soil microclimate. For instance, it has been observed that the abundance and diversity of woodlice, snails, and beetles increases toward the lower parts of dolines, while the surrounding plateau possesses higher species richness and abundance of ants and spiders in forested dolines (Sólymos et al., 2009; Vilisics et al., 2011; Kemencei et al., 2014; Bátori et al., 2022). In the present study, except for open doline (2), Collembola showed overall increasing values of community parameters toward the bottoms of the dolines, which were characterized by generally low soil temperatures and low carbon and nitrogen content. This pattern was supported by significantly negative correlations between the Shannon diversity index and minimum soil temperatures as well as carbon content. Moreover, dolines with this community pattern were characterized by a pronounced gradient of vegetation, with their bottoms occupied by predominantly hygrophilous associations.

The observed pattern of the soil Collembola communities differed considerably from that of a collapse doline with a much steeper microclimatic gradient, the entrance slope of the Silická ľadnica Ice Cave. There, a decreasing trend in species richness and an increasing trend in mean abundance of the Collembola were found toward the cold part of the gradient, resembling a typical latitudinal pattern of abundance and species richness (Raschmanová et al., 2013; Potapov et al., 2023). This indicates that the geomorphology of the karst dolines considerably shapes the distribution and structure of soil fauna communities.

Forest age, affecting the quantity of light, the microclimate, or soil nutrients, may also shape/influence the species richness and community composition of soil arthropods; in addition, this can be a good indicator for the number of cold-adapted species in dolines (Bátori et al., 2023a). The vegetation effect may be direct by influencing the organic matter layer or indirect by changing the soil microclimate and quality or quantity of plant litter (Eaton et al., 2004). In this study, the forested dolines, characterized by the presence of an approximately 70-year-old oak forest, showed markedly higher diversity and abundance of Collembola and hosted higher numbers of cold-tolerant species compared to the open dolines, which can be attributed to (1) the presence of leaf litter offering Collembola a high microhabitat heterogeneity and an indirect food source (Hopkin, 1997; Rusek, 2001; Bardgett, 2002; Salamon et al., 2004; Perez et al., 2013) and to (2) more favorable microclimatic conditions (i.e., lower soil temperatures and temperature fluctuations).

Besides the microclimate, the content of soil organic carbon is also an important factor, determining the diversity and composition of the Collembola communities (Rendoš et al., 2016; Jureková et al., 2021). However, in this study, soil nutrients had a less evident effect on community parameters than soil microclimate, showing the higher values of community parameters at sites associated with the low content of soil organic carbon.

Two types of karst dolines hosted different communities of soil Collembola. Within the open dolines, a clear delimitation concerning the microclimate was observed at warm sites associated with the xerothermophilous species Isotoma anglicana and Proisotomodes bibunctatus and the thermophilous and mesophilous Protaphorura serbica. In contrast, colder sites were associated with the eurythermic and mesophilous Lepidocyrtus lignorum, Parisotoma notabilis, and Megalothorax minimus. The forested dolines did not show such a visible trend, which reflects the lower soil temperature fluctuations in this environment.

Implementing functional trait approaches in ecology is a suitable way of understanding the response of soil arthropod communities to environmental gradients, as was shown in recent studies by Bátori et al. (2019a, 2022) and Marcin et al. (2022). These studies documented the potential of karst dolines to facilitate the persistence of distinctive arthropod functional/ecological groups on a relatively small scale. It was observed that open karst dolines can provide microhabitats for various functional groups of soil Collembola, where hygrophilous and mesophilous species preferred the bottoms of open dolines, and xerophilous/xeroresistant and thermophilous species clearly preferred the warm sites of the dolines, that is, the S-facing slopes and the plateau sites (Marcin et al., 2022). A significant preference of thermophilous and xerophilous/xeroresistant ecological groups for warm sites was observed in both the open and forested karst dolines. Moreover, more marked distributional patterns of cold-adapted and hygrophilous species were observed within the forested dolines, as they strongly preferred cold N-facing slopes and wet doline bottoms, respectively. This strong pattern of preference of ecological groups for microclimatically different sites along these dolines points out that despite the buffered climatic conditions typical for forest habitats, topographic and microclimatic conditions can strongly delimit the distribution of ecological groups of soil Collembola not only in open but also in forested dolines. Topography can indeed strongly affect the functional composition of Collembola communities, even if without a visible effect on species composition and diversity (Hishi et al., 2022).

At present, the threat of global warming is placing numerous species at risk of extinction due to shifts in local temperature and rainfall patterns and thus unsuitable climatic conditions for many species (Suggitt et al., 2018). However, topographically diverse habitats can modify or mitigate these effects on a regional level by generating distinct meso- or microclimatic gradients. These unique climate conditions can form a complex assortment of suitable habitats, allowing many species to endure unfavorable periods of climate conditions (Raschmanová et al., 2008; Bátori et al., 2014; Raschmanová et al., 2018; Suggitt et al., 2018; Papáč et al., 2019; Marcin et al., 2021; Bátori et al., 2022; Marcin et al., 2022). In this study, a considerably high number of endemic species were found especially in the forested dolines, with most of them being cold-adapted species. Our study pointed out that karst dolines serve as natural islands of the exceptional diversity of soil fauna and provide important microrefugia for endemic species (Raschmanová et al., 2008; Vilisics et al., 2011; Bátori et al., 2014; Raschmanová et al., 2018; Bátori et al., 2019b). The presence of endemic species and species characterized by narrow distribution ranges seriously outline the importance of karst dolines from a conservation point of view. Numerous studies have shown that the cold environments of some karst landforms maintain high microclimatic stability (Maclean et al., 2015; Yao et al., 2016; Kováč, 2018; Perşoiu and Lauritzen, 2018; Suggitt et al., 2018; Liu et al., 2019). The high number of cold-adapted species and their marked distributional patterns in forested dolines highlight the favorable microclimatic conditions for such specialized cold-tolerant species, more specifically cold N-facing slopes and doline bottoms.

Both types of studied dolines could potentially provide microclimatic refugia for cold-adapted species by providing relatively cool and humid conditions at the bottoms and poleward-facing slopes (Bárány-Kevei et al., 2021), which is more pronounced in forested dolines due to the presence of forest cover, resulting in a buffering of the local microclimate against the surrounding climate (Wallwork, 1970; Coleman and Crosley, 1995; De Frenne et al., 2019; De Frenne et al., 2021). The gradient of microclimatic conditions resulting from karst landscape topography and environmental stability could support the long-term persistence of at least some species, such as cold-adapted/psychrophilous species. Therefore, the assessment of the function of the dolines as potential climate refugia for psychrophilous forms of soil fauna is crucial for improving current and developing future conservation strategies, especially in the context of global climate change (Dobrowski, 2011; Raschmanová et al., 2018; Bátori et al., 2019a; Aguilon et al., 2020; Kiss et al., 2020).

Our findings highlight enclosed surficial karst dolines as landforms of great microclimatic and habitat heterogeneity. The forested dolines were characterized by a markedly buffered microclimate with lower differences in soil temperature between topographically distinct sites and provided conditions for markedly higher species richness and abundance, with increasing patterns toward the inner parts of the dolines. The forested dolines hosted high numbers of vulnerable taxa, such as cold-adapted species, thus documenting that karst landforms serve as local refugia of specialized (cold-tolerant) species.

Both types of dolines showed the capability to maintain distinctive ecological groups of Collembola reflecting microclimatic partitioning at sites. More marked distributional patterns of the cold-adapted and hygrophilous species within the forested dolines were observed. More than one-third of Collembola species occurred only in the inner parts of both types of dolines, pointing out their role in facilitating the persistence of some invertebrates, thus importantly supporting local biodiversity.

Considering the special topographic and microclimatic conditions, soil fauna patterns, and species composition, forested and open karst dolines are especially important and valuable for nature conservation.

All data generated during the study are presented in this article.

The authors thank colleague Ľuptáčik (Pavol Jozef Šafárik University in Košice, Košice, Slovakia) for assistance during the fieldwork. R. Šuvada (Administration of the Slovak Karst National Park, Brzotín, Slovakia) is acknowledged for the analysis of the vegetation associations of the studied dolines. They are grateful to D. L. McLean for the linguistic correction of the manuscript.

This work was funded by the Slovak Research and Development Agency (grant APVV-21-0379) and by the Slovak Scientific Grant Agency (grant VEGA 1/0438/22).

The authors declare no competing interests.

Contributed to conception and design: MM, NR, ĽK.

Contributed to acquisition of data: MM, NR.

Species identification: MM, NR.

Contributed to analysis and interpretation of data: DM, JŠ, JK.

Drafted and/or revised the article: MM, NR, DM, ĽK.

Approved the submitted version for publication: MM, NR, DM, JŠ, JK, ĽK.