Karst dolines provide diverse microhabitats for different functional groups in multiple phyla

Fine-scale topographic complexity creates important microclimates that can facilitate species to grow outside their main distributional range and increase biodiversity locally. Enclosed depressions in karst landscapes (‘dolines’) are topographically complex environments which produce microclimates that are drier and warmer (equator-facing slopes) and cooler and moister (pole-facing slopes and depression bottoms) than the surrounding climate. We show that the distribution patterns of functional groups for organisms in two different phyla, Arthropoda (ants) and Tracheophyta (vascular plants), mirror this variation of microclimate. We found that north-facing slopes and bottoms of solution dolines in northern Hungary provided key habitats for ant and plant species associated with cooler and/or moister conditions. Contrarily, south-facing slopes of dolines provided key habitats for species associated with warmer and/or drier conditions. Species occurring on the surrounding plateau were associated with intermediate conditions. We conclude that karst dolines provide a diversity of microclimatic habitats that may facilitate the persistence of taxa with diverse environmental preferences, indicating these dolines to be potential safe havens for multiple phyla under local and global climate oscillations.


Results
Microclimate. Temperatures were higher on south-facing slopes (T 24 = 25.9 °C; T d = 33.8 °C; T n = 13.8 °C) than in other microhabitats (Figs 1 and 2B). The mean daily temperature was similar on north-facing slopes (T 24 = 20.1 °C) and in doline bottoms (T 24 = 20.6 °C), while night temperatures were lowest in bottoms (T n = 9.78 °C). Night temperatures were similar on north-facing slopes (T n = 11.5 °C) and the plateau (T n = 11.8 °C). The mean daytime temperature was higher in bottoms (T d = 27.7 °C) than on north-facing slopes (T d = 25.7 °C). Temperatures on the plateau were intermediate (T 24 = 24.0 °C; T d = 31.8 °C; T n = 11.8 °C). Mean daily relative humidity was lowest on south-facing (RH 24 = 68.0%) and highest on north-facing slopes (RH 24 = 78.4%). The values were intermediate in bottoms (RH 24 = 76.2%) and on the plateau (RH 24 = 70.0%). At night, relative humidity was similar (RH n = 92.9-94.0%) in all microhabitats. However, daytime relative humidity was usually higher on north-facing slopes (RH d = 68.2%) than in bottoms (RH d = 65.2%), on south-facing slopes (RH d = 51.7%) and on the plateau (RH d = 55.1%). A small, intermittent temperature increase in all microhabitats was recorded around 1:30 hours and indicated an inflow of warmer air from the surrounding lower altitudes. species composition. A total of 14 ant (nine from baiting and another five from hand collecting) and 145 plant species (from the plots) were recorded in our sites ( Fig. 1G; Supplementary Tables S1 and S2). Hand collection yielded 52 ant individuals, 13 were found on south-facing slopes, five in bottoms, 11 on north-facing slopes and 23 on the plateau. In terms of diagnostic species, south-facing slopes had two ant and 13 plant species, while two ant and 20 plant species were identified for the plateau. North-facing slopes and bottoms had 15 and nine plant species, respectively, and one ant species ( Fig. 2A; Table 1). Dry grasslands dominated the south-facing slope of dolines 1 and 3, while semi-dry grasslands covered the south-facing slope of doline 2 and the major part of the plateau. North-facing slopes and doline bottoms were dominated by wet meadows.
NMDS ordinations of bait stations (stress factor: 0.119) and vegetation plots (stress factor: 0.064) showed that the species composition of ant and plant assemblages differed among microhabitats (Fig. 3A,B). These differences were significant, except between ant assemblages of bottoms and north-facing slopes (PERMANOVA; Table 2). Plant assemblages of all microhabitats displayed highly significant differences (p < 0.001). T 24 , T d and RH d were significantly related to the ordination of ant, and T 24 , T d , RH 24 and RH d to the ordination of plant assemblages (Supplementary Table S3), with assemblages on south-facing slopes being associated with higher temperatures and lower humidity, and assemblages on north-facing slopes and in doline bottoms being related to lower temperatures and higher humidity.  Tables S4 and S5). On north-facing slopes, 64% of the hand-collected ants (individuals) were species adapted to cooler conditions, and 82% of the individuals were species adapted to intermediate moisture conditions. In other microhabitats, most of the hand-collected ants (100% on south-facing slopes, 80% in bottoms and 52% on the plateau) were species adapted to warmer and drier conditions. Ant (collected from bait stations) and plant species adapted to warmer and/or drier conditions occurred more frequently on south-facing slopes than on north-facing slopes and in doline bottoms (Figs 4 and 5;  Supplementary Tables S4 and S5). Conversely, ant and plant species adapted to cooler and/or moister conditions were generally most frequent on north-facing slopes and in doline bottoms compared to south-facing slopes and the plateau. We did not find any ant species adapted to moister conditions and no plant species adapted to cooler conditions on south-facing slopes.

Discussion
Karst dolines are topographically complex environments that provide a variety of microhabitats. We have demonstrated that these habitats can be cooler and moister or warmer and drier than the surrounding plateau, providing high environmental heterogeneity at very fine scales. We further documented that the fine-scale distributions of functional groups in two different phyla (Arthropoda and Tracheophyta) correspond to these different microhabitats. To our knowledge, this is the first study to illustrate that the fine-scale topography of dolines provides microhabitats for diverse functional groups (cool-and moist-adapted versus warm-and dry-adapted) of both ants and vascular plants within tens of meters. www.nature.com/scientificreports www.nature.com/scientificreports/ Topographic complexity increases the climatic variability in an area over fine scales 24,51,52 . Despite recording temperature and relative humidity for only a relatively short period of time (24 hours), we found that karst dolines introduce great variation in microclimates. For instance, mean daytime temperatures were more than 8 °C warmer on south-than on north-facing slopes (Fig. 2B). Previous microclimatic studies over longer time periods, i.e. from a few days to a year, also indicate north-facing slopes and bottoms of dolines to be cooler and more humid than south-facing slopes and the surrounding plateau 23,35,[53][54][55][56] . In addition to receiving less solar radiation, doline bottoms tend to receive more water 57 , likely contributing to the higher relative humidity and lower temperatures recorded in these habitats 58,59 . North-facing slopes and bottoms of dolines may retain snow cover longer than south-facing slopes and plateaus 60 .
Fine-scale variation in environmental conditions affects the distributions of organisms [61][62][63] . Both ant and vascular plant species responded to microclimatic variation in our study. The cooler and moister north-facing slopes and bottoms of dolines in Bükk acted as key habitats for plants adapted to cooler and/or moister conditions (e.g. Bupleurum longifolium and Iris sibirica) and ants adapted to cooler conditions (e.g. Myrmica ruginodis) ( Fig. 2A). On the other hand, south-facing slopes acted as key habitats both for ants (e.g. Myrmica sabuleti) and plants www.nature.com/scientificreports www.nature.com/scientificreports/ (e.g. Iris variegata) adapted to warmer and/or drier conditions, while many ant (e.g. Lasius niger) and plant species (e.g. Galium verum) found on the plateau indicated intermediate temperature and/or moisture conditions ( Fig. 2A). The ant species recorded are not known to form strict trophic relationships with plants, and the main drivers of the observed patterns seem to be their temperature and moisture preferences. Our results underline the importance of considering fine-scale environmental variation when investigating the distribution of biodiversity 63,64 . Even the highest resolution (50-100 m) climate data currently available 7,13 would likely be insufficient to detect the environmental heterogeneity provided by the karst dolines in our study, which are <200 m wide and <20 m deep.
We have documented a strong, concerted response to fine-scale topography by the distributions of microclimate and species in two major taxonomic groups (ants, Arthropoda; and vascular plants, Tracheophyta). However, other taxa also display distribution patterns reflecting changes in environmental conditions in dolines over short distances. For instance, the increased soil moisture content on north-facing slopes and bottoms of solution dolines in the Aggtelek Karst area (Hungary) has been shown to provide suitable habitats for several cool-adapted species of land snails (Mollusca) 49 and terrestrial isopods (Arthropoda) 65 . Dolines in Mexico, Australia and China ('cenotes' and 'tiankengs') have also been shown to maintain populations of rare taxa in various phyla 46,66-68 , highlighting dolines as important safe havens for biodiversity on a global scale.
Species may respond to climate changes by range-shifting 69 or by persisting in environmentally stable habitats 70,71 . In karst dolines, species could potentially track suitable microclimates over a wide range of regional climatic changes with minimal movement because cooler/moister and warmer/drier microclimates vary over very short distances. As few areas can buffer opposing trends in environmental conditions 72 , karst dolines may be particularly important for maintaining biodiversity through time. Therefore, they could be considered high-capacity microrefugia 73 . The highest-capacity microrefugia for cool-adapted taxa can usually be found in cold, humid and topographically complex areas 74,75 .
Therefore, karst dolines may play important roles in facilitating the persistence of different phyla under global warming, which poses a serious threat to global biodiversity 76 . Regional predictions of climate change suggest that warming in East-Central Europe will continue in the coming decades 77 . These changes are already impacting the distributions of ants 78 and vascular plants 79 in a sand-dune area in Hungary, with drought-tolerant species replacing dune slack species over the last decades. Therefore, species adapted to warmer and/or drier conditions are expected to expand their distribution from south-facing slopes of dolines to surrounding areas. However, north-facing slopes and bottoms of dolines could provide important microrefugia from global warming by facilitating the persistence of species adapted to cooler and/or moister and to intermediate conditions. The retention of cooler microclimates in these habitats may be facilitated by lower solar radiation, thicker soil layer, higher soil moisture and cool-air pooling 27,58,80 .
Although our data supports dolines to be safe havens for relict plant species (e.g. Aconitum variegatum, Bupleurum longifolium and Dracocephalum ruyschiana in Bükk) in the current climate, future studies should aim to confirm the status of dolines as refugia (i.e. places providing environmental conditions that are comparatively stable over long time periods) under global warming. While our 24-hour data demonstrate that dolines are currently providing cooler and warmer microclimates than the surrounding plateau, this does not necessarily prove stability and long-term monitoring would be needed for this. Available data from north-facing slopes in the Northern Hemisphere does suggest that such habitats undergo slower changes under global warming 81,82 . Alternatively, the microclimate of dolines could be investigated along a temperature gradient using a space-for-time substitution approach to determine if north-facing slopes and bottoms of dolines indeed retain  www.nature.com/scientificreports www.nature.com/scientificreports/ more stable microclimates. Finally, functional traits can be reflective of long-term environmental stability and therefore could provide important eco-evolutionary information about refugia 83 .
We conclude that enclosed depressions in karst surfaces provide a diversity of microclimates that have the potential to enable the persistence of various taxa in different phyla and under various climatic trends. These dolines may be vital for facilitating the in situ persistence of numerous species under local and global climate oscillations. This implies that modelling of climate change impacts on the distribution of biodiversity will need to consider fine-scale topographic variation occurring within tens of meters to arrive at accurate predictions.  www.nature.com/scientificreports www.nature.com/scientificreports/ Methods study area. Our study area was located on the karst plateau of the Bükk Mountains (48°04′31″N, 20°29′57″E), in northern Hungary, at an altitude of approximately 780 m (Fig. 1A,B) in the beech (Fagus sylvatica) forest zone. This mountain range is believed to be an important refugial area in Hungary, supporting relict plant populations from both warmer (e.g. Clinopodium thymifolium, Cotinus coggygria and Ferula sadleriana) and cooler (e.g. Aconitum variegatum, Bupleurum longifolium and Dracocephalum ruyschiana) periods 84,85 . The plateau has a cool and humid climate, with a mean annual temperature of 6.3 °C and a mean annual precipitation of 800 mm. The plateau has typical karst landform features, such as solution dolines [86][87][88] , with a bowl-shaped geometry (Fig. 1C) and unique microclimate. At night, cold-air pooling occurs in these depressions, and the occurrence of frost www.nature.com/scientificreports www.nature.com/scientificreports/ or fog is possible all year round 55 (Fig. 1D). According to some researchers 89 , the coldest areas in Hungary can be found in the non-forested dolines of Bükk. Previous investigations showed that north-facing slopes of these dolines are consistently cooler and moister than the surrounding plateau, and that south-facing slopes provide the warmest microhabitats 54,90 . The study area is well known for its unique wildlife and is part of the strictly protected area network of the Bükk National Park. The entire area is fenced to prevent overgrazing, soil erosion and the illegal collection of wild plants. Semi-dry grasslands and wet meadows are the dominant vegetation types within the fenced area 91 (Fig. 1E,F). www.nature.com/scientificreports www.nature.com/scientificreports/ sampling design. Three large solution dolines were selected (Fig. 1B,C). Dolines 1 and 2 had diameters of 100 and 70 m and depths of 17 and 15 m, respectively. The longer diameter of doline 3 was 190 m, while the shorter was 65 m, its depth was 13 m. Six sampling sites were established per each doline (18 sites in total), one on the south-facing slope, one in the bottom, one on the north-facing slope and three on the surrounding plateau. The sites were located at least 20 m from each other (Fig. 1G). Sampling and microclimate measurements were carried out in August, under clear weather conditions, at the peak of the growing season.
Ants and plants were selected as focal taxa because ant colonies and plants share many similarities 92 . Both groups usually 'nest' in or on the ground and use their modules (e.g. plant roots and ant workers) to forage in the surrounding habitat 93,94 . In addition, due to the relatively fixed location of ant colonies and plants, competition in both groups is confined to well-defined zones. Similarities also exist in their functional roles in a given community (e.g. subordinate, specialist and cryptic species). Finally, ant foundresses (i.e. colony-founding queens) can be considered analogous to dispersing plant seeds 93,95 .
Considering that the study area is a strictly protected nature reserve, we used only non-destructive sampling methods such as baiting and hand collecting to assess the species diversity and relative abundance of ants. At each site, we placed five bait stations in a cross-shaped pattern at 5-m intervals (90 bait stations in total) (Fig. 1H). Baits were plastic discs (8 cm in diameter) with a quarter-teaspoon of a mixture of tuna and honey as a food reward. Foraging activity on baits was monitored every 40 minutes from 7:00 to 9:40, overlapping with the daily period of peak ant activity. During each observation, we recorded the presence and number of workers of each ant species on the bait. Baits were replenished as necessary. In addition, we also performed hand collecting to sample those ant species that may have not visited the baits. We visually searched the ground surface in each site for 5 minutes, hand collecting any individuals (workers, incipient queens, etc.) found. Ants were identified to morphospecies or genus level in the field, and representatives were collected and preserved in 95% ethanol for later species determination 96 . Field-collected specimens were identified in the laboratory using the keys of Seifert 97 and Czechowski et al. 98 . All the collected specimens were deposited at the Department of Ecology, University of Szeged. We used Bolton's catalogue 99 and the Hymenoptera Name Server 100 to determine the valid names of all ant species.
For plants, five randomly selected 1 m × 1 m plots were established in each site (90 plots in total). We recorded the presence/absence of all vascular plant species in all plots. Nomenclature follows 'The Plant List' (www.theplantlist.org).
To provide information on the microclimate of the study area, air temperature (T) and relative air humidity (RH) were recorded every minute for 24 hours using Voltcraft DL-121TH data loggers. Sensors were suspended 10 cm above the ground to allow sufficient wind flow to ensure that no humidity was trapped by the sensor casing and actual air temperature and humidity were measured. species grouping. We classified all ant and plant species according to their temperature and moisture requirements following the methods of Czechowski et al. 98 Table S6 and Table S7). These were reduced to six functional groups applied to both ants and plants: (1) Supplementary Tables S6 and S7). We did not analyse combined groups (temperature + moisture) because that would have made the interpretation of results difficult (with many groups), especially in the case of plants. All six main groups of plants and five main groups of ants (groups 1-4 and 6; ants adapted to moister conditions were absent in our study) were analysed. Because the temperature and moisture requirements of ant and plant species are not the same, direct comparisons between them were not possible. Data analysis. The temperature and relative humidity data were averaged over 10-minute intervals across all sites of south-facing slopes, bottoms, north-facing slopes and the plateau, respectively, and plotted using a line graph. Extreme environmental values are generally more informative regarding the distribution of organisms, but maximum relative humidity values in our case often reached 100%, therefore we considered mean values more suitable for differentiating between the microclimatic properties of microhabitats than extreme values.

for ants and Borhidi 101 for plants (Supplementary
From the site-averaged data, we calculated the mean daily temperature (T 24 ) and relative humidity (RH 24 ), mean daytime temperature (T d ) and relative humidity (RH d ), and mean night temperature (T n ) and relative humidity (RH n ). We also calculated these microclimate variables separately for each site, and used them in multivariate analyses.
The diagnostic ant and plant species of the microhabitats were determined by calculating the phi (Φ) coefficient of association between species and habitat. Calculations for ants were based on data from bait stations. Species with Φ > 0.2 were considered diagnostic for ants and species with Φ > 0.3 for plants. Different threshold values for ants and plants were used due to the differences in the total number of species within each of these taxonomic groups. Non-diagnostic species were excluded with Fisher's exact test (p < 0.05) following Tichý and Chytrý 102 . Fidelity measures were calculated using the JUICE program 103 .
We used permutational multivariate analysis of variance (PERMANOVA) to test the effect of microhabitats (south-facing slopes, bottoms, north-facing slopes and the plateau) on the species composition of ant and plant assemblages. We used the raw presence/absence data of species for each sampling plot in the source matrices. We applied the Jaccard dissimilarity, performed 5000 permutations and also accounted for the nested design of the data set. When a PERMANOVA yielded significant results, we calculated pairwise PERMANOVAs among the microhabitat types. PERMANOVAs were calculated in R 104 using the adonis function of the 'vegan' package 105 . We used the FDR (false discovery rate) method to adjust p values for multiple comparisons (p.adjust function). We prepared non-metric multidimensional scaling (NMDS) ordinations to visually illustrate compositional differences. To remove the confounding effect of the nested data structure on the resulting point clouds, we lumped Scientific RepoRts | (2019) 9:7176 | https://doi.org/10.1038/s41598-019-43603-x www.nature.com/scientificreports www.nature.com/scientificreports/ data from the different sampling plots of each site of each doline together and used the frequency data, ranging from one to five, of the species in the source matrices. We used Euclidean distances and two dimensions (after assessing stress factors for one to five dimensions). NMDS ordinations were done using the metaMDS function of the 'vegan' package. To assess the relationships between microclimate variables (T 24 , RH 24 , T d , RH d , T n and RH n ) and species assemblages, we fitted environmental vectors onto the ordination space using the envfit function and calculated correlations between ordination values and fitted vectors.
We used mean-and-whisker plots to illustrate the distribution of the various functional groups in different microhabitats. To test if these differences were significant we used generalized linear mixed-effects models (GLMM). Calculations for ants were based on data from bait stations. Five models were built for ants (three for temperature and two for moisture) and six models for plants (three for temperature and three for moisture). In the full models, different microhabitats were included as fixed factors, the number of ant and plant species as dependent variables, and location (i.e. doline 1, 2 and 3) as the random factor. We transformed the data of ants to binary scale (presence/absence) and used a binomial error term because each functional group had a high preference for one or a few sites and were very rare in other sites, leading to zero inflation of the data. No transformation was needed for plants, and we used Poisson or, if overdispersion was detected, negative binomial error term. GLMMs were performed in R using the glmer function of the 'lme4' package 106 . Full models were tested for significance with analysis of variance, using the Anova function of the 'car' package 107 . Pairwise comparisons of factor levels were undertaken with the relevel function and the FDR method (p.adjust function) was used to correct p values for multiple comparisons.

Data Availability
The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.