Introduction

Glaciers and ice sheets cover ~ 10% of the global surface and have attracted geomorphologists and climatologists as one of the significant drivers of Earth's landscape and climate (Lewandowski et al. 2020; Schuler et al. 2020). Although historically considered unfriendly for life, recent biological research provides us with evidence that glaciers also constitute important component of a global ecosystems. They are inhabited by a myriad of primarily microscopic organisms forming relatively short food chains comprising single celled as well as multicellular organisms (Hodson et al. 2008; Anesio and Laybourn-Parry 2012; Cook et al. 2016; Zawierucha et al. 2021). The most biodiversity-rich habitats on glaciers are cryoconite holes (Fig. 1), water-filled reservoirs with dark sediment at the bottom that form a suitable habitat for bacteria, protozoans, plants, and animals (Takeuchi et al. 2010; Zawierucha et al. 2021; Rozwalak et al. 2022).

Fig. 1
figure 1

Ebbabreen and cryoconite holes: A Ablation zone of Ebbabreen; B Surface of the ablation zone on Ebbabreen; and C, D Cryoconite holes on Ebbabreen

Full size image

Many organisms inhabiting glaciers and ice sheets represent specialized lineages, which often possess morphological, physiological, and behavioral adaptations to survive permanently low temperatures, strong seasonality, frequent freeze–thaw cycles, and high doses of UV radiation (Dastych et al. 2003; Singh et al. 2014; Dial et al. 2016a; Shain et al. 2016; Zawierucha et al. 2019a). A question, therefore, emerges as to whether their demes form genetically interconnected metapopulations spread across separated glaciers or retain locally restricted genetic variants potentially representing local adaptations (Zawierucha et al. 2023).

Given their high dispersal potential, it is currently debated whether the microscopic animals follow the same biogeographical rules of unicellular microorganisms where some taxa have a really global range (e.g., red snow algae or cyanobacteria; Segawa et al. 2017, 2018) or have a number of range-restricted species as observed for the macrofauna (e.g., Cho and Tiedje 2000; Whitaker et al. 2003; Darling et al. 2004; Foissner 2005; Telford et al. 2006; Fontaneto and Hortal 2008; Fontaneto 2011). Biogeography of cryophilic organisms is supposedly further shaped by special spatio-temporal characteristics of glacier habitat. On a large scale, especially northern glaciers represent geographically fragmented refugia hosting populations of cryophiles that once might have been part of much larger ancient populations during glacial periods. On a small scale, each isolated glacier represents a highly heterogeneous and dynamic structure where individual cryoconite holes may remain isolated even on multi-seasonal scales, like in the Antarctic glaciers (Darcy et al. 2018), but can also represent highly dynamic structures that are repeatedly flushed and likely recolonized over several weeks or a single day (Takeuchi et al. 2018; Zawierucha et al. 2019b). Glacial biogeography also has an important temporal aspect since organisms that inhabit the terrestrial cryosphere may become “conserved” in permanently frozen ice layers and, owing to their ability to form cryobiotic stages, be viable for decennia or longer (Shmakova et al. 2021).

While biogeography of glacial metazoans received less attention than primary producers and heterotrophic bacteria, several studies revealed interesting patterns, ranging from allopatric fragmentation to long-distance dispersals across several hundreds of kilometers (e.g., Dial et al. 2016b; Shain et al. 2016; Zawierucha et al. 2018b; Hotaling et al. 2019).

A prominent metazoan group found on glaciers is the phylum Tardigrada (also known as water bears), whose cryptobiotic state enables them to survive harsh conditions (Persson et al. 2011; Jönsson et al. 2016). Water bears are common in polar ecosystems and dominate on glaciers, where they sometimes represent the only group of metazoans (Convey and McInnes 2005; Zawierucha et al. 2021). Tardigrades recorded from cryoconite holes appear sensitive to high temperatures, but they are able to withstand freezing or desiccation at low temperatures (Zawierucha et al. 2018b, 2019b, a).

Recent taxonomic work led to the description of new species from glacial habitats that range from tropical regions (Zawierucha et al. 2018a) to high latitudes, e.g., the Pilatobius glacialis Zawierucha, Buda & Gąsiorek, 2020, which dominates cryoconite holes across Svalbard Archipelago (Zawierucha et al. 2020), or the cryophilic specialist genus Cryoconicus (Ramazzottiidae: Eutardigrada) from Central Asian glaciers (Zawierucha et al. 2018b). Cryoconicus possesses a suite of interesting adaptations to conditions persisting on glaciers, i.e., ability to withstand low temperatures, including being frozen for 11 years, or dark pigmentation that is hypothesized to protect it from intense UV irradiation (Zawierucha et al. 2018b). Especially Cryoconicus kaczmareki Zawierucha et al. 2018b appears to be a true glacier specialist reported from several Central Asian glaciers (Zawierucha et al. 2018b). Analysis of its mtDNA variability revealed that (i) identical COI haplotypes were shared among some Asian populations separated by ca. 1000 km and (ii) even diverged haplotypes were found within the same Asian populations. These patterns imply either good dispersal abilities of the species or recent fragmentation from a previously larger population (Zawierucha et al. 2018b, 2019a).

To investigate the biogeography of Cryoconicus and to test its potential presence on glaciers worldwide, we examined our extensive collection of cryoconite material and studied the literature on tardigrades from glaciers across the world (Fig. 2). We applied morphological and molecular analysis of newly found specimens in the Arctic to confirm their conspecificity with described Asian and Antarctic Cryoconicus species. Further, we employed phylogeographic methods to investigate population connectivity and dispersal history across the distribution ranges. As Cryoconicus-like taxa were newly found on restricted glacial outlets of the Spitsbergen ice cap, we included additional seasonal resampling of some of the glaciers to test the multiyear stability of local populations. Finally, we compared the population genetic patterns and local abundances of Cryoconicus species with those of Pilatobius glacialis, a dominant tardigrade on Svalbard glaciers, in order to test whether the observed abundance ranking of glacial species remained stable in the long term or underwent fluctuations in the recent past.

Fig. 2
figure 2

Map of the World showing the distribution of C. kaczmareki and C. antiarktos and sampled glaciers (empty ellipses and dots) and other habitats (green circles). Numbers correspond to sites as mentioned in Online Resource 1. Black & empty ellipses and circles—sampled glaciers without Cryoconicus species; blue & empty ellipses and circles—terra typica of C. kaczmareki; blue circles with yellow asterisk—new records of Cryoconicus in Svalbard and Asia; brown circles—presence of C. antiarktos; full green with asterisk—bryophyte with a new record of C. antiarktos; full green—original locus typicus of C. antiarktos. Details are provided in Online Resource 1

Full size image

Materials and methods

Comparison to worldwide data collection

To screen for the presence of tardigrade species belonging to the genus Cryoconicus on glaciers worldwide, we analyzed available literature data focusing explicitly on tardigrades from 60 glaciers with differing characteristics: (i) geological settings, (ii) climate, (iii) type of glacier (e.g., tidewater, piedmont, valley, ice sheet), (iv) glacier thermal regime (i.e., polythermal, cold-based, temperate), (v) light availability (e.g., seasonal cycles in high latitudes, daily cycles in lower latitudes), and (vi) elevation (e.g., marine-terminating tidewater glaciers, ice caps, and high mountain glaciers up to 5895 m a.s.l.). In some studies, cryoconite habitats were sampled during more than one sampling campaign over multiple years (e.g., Ebbabreen, Longyearbreen, Russell Glacier, Blåisen Glacier, Forni Glacier, Ecology Glacier, Canada Glacier), while other glaciers were sampled only once (e.g., Kersten Glacier, La Conejeras Glacier, Marr Glacier). A list of sampled glaciers is provided in Online Resource 1. The presence/absence of Cryoconicus was evaluated using literature and original data (utilizing both classical [microscopic observation] and molecular approaches). In these studies, cryoconite was sampled either using independent sterile disposable Pasteur pipettes, or by scoops, sterilized before taking each sample. Samples were placed directly in vials, jars, or sterile Whirl–Pak® bags. All sediment samples were either immediately frozen after fieldwork (a few hours), frozen after a few days (samples from Baltoro Glacier), or preserved in 70% or 96% ethyl alcohol for transport to laboratories.

Additionally, we analyzed the morphology and DNA of species of Cryoconicus found in a moss sample from the Ablation Valley in Alexander Island, Antarctica, in order to increase the dataset of limno-terrestrial Cryoconicus species and knowledge on its distribution.

Tardigrade collection and DNA barcodes from individual specimens

Tardigrades were extracted from cryoconite samples under stereomicroscopes. A limited number of specimens belonging to the genus Cryoconicus were recovered from sampling on two Arctic glaciers (Ebbabreen and Nordenskiöldbreen, see the Results). Some of these animals were preserved for detailed microscopic analysis and the rest was used for DNA analysis. Genomic DNA was successfully extracted from thirteen individuals from Ebbabreen and a single individual from Nordenskiöldbreen using the DNeasy Blood and Tissue Kit (Qiagen GmbH, Hilden, Germany) according to the manufacturer’s protocol. To retrieve the exoskeletons after digestion with Proteinase K and lysis, only the supernatant was transferred to the spin columns with a silica-based membrane. Retrieved exoskeletons were subsequently mounted in Hoyer medium for morphological determination. One mitochondrial (cytochrome-c oxidase subunit I, COI) and one nuclear (D1-D3 region of 28S rDNA) gene fragments were amplified and sequenced using protocols and primer described in Online Resource 2.

To compare the new samples of Cryoconicus with their morphologically determined conspecifics, we used previously published COI and 28S rDNA sequences of C. kaczmareki (sequences No. MG432803-9 and MG432797.1) and Cryoconicus antiarktos Guidetti et al. 2019 (MK879525-31, MK879544.1 and MK879543.1). Sequences were visually inspected in BioEdit (Hall and Hall 1999) for accuracy of base calls, presence of potential heterozygotes and stop codons. Alignment was performed using the ClustalW program with the default settings for the COI gene and using MAFFT for the 28S rDNA dataset (Nakamura et al. 2018).

For comparative purposes, we also used the same protocols to extract DNA and amplify the COI barcode of P. glacialis, the dominant species on Arctic glaciers, which was collected from four glaciers draining to Isfjorden, Central Svalbard (Terra typica of species), i.e., the same biogeographical zone for the Svalbard population of C. kaczmareki (Table 1).

Table 1 Diversity statistics of Cryoconicus kaczmareki and Pilatobius glacialis
Full size table

Analysis of intraspecific diversity and molecular dating

Summary statistics, including haplotype, nuclear diversities, and their variance, were calculated in DNAsp v. 5 (Librado and Rozas 2009). Fine-scale phylogenetic resolution within and among C. kaczmareki and P. glacialis populations was obtained with the parsimony tree reconstruction using the TCS algorithm (Clement et al. 2000) as implemented in PopArt (Leigh and Bryant 2015).

To estimate the divergence of the Arctic population of C. kaczmareki from the Asian population, we used coalescent analysis implemented in BEAST 1.8.0 (Drummond et al. 2012), running the analyses on both, all codon positions, and on 3rd codon positions only, in order to minimize the effects of purifying selection (e.g., Černá et al. 2017). Models of sequence substitution were estimated in W-IQ-Tree (Trifinopoulos et al. 2016) and set to TN with either invariant sites or with empirical base frequencies, respectively. Trees were constructed using a coalescence model of constant size under the strict molecular clock assumption and were estimated with two independent runs of 1 × 108 iterations, sampling every 10,000th generation. The results were checked for convergence in Tracer v1.5 (Rambaut et al. 2009) and combined in the LogCombiner program using 10% burn-in.

Since the tardigrade mtDNA mutation rate is unknown, we applied two different approaches to translate the node depths into absolute times. First, the COI tree based on 3rd positions was calibrated with the neutral mutation rates estimated for Drosophila melanogaster Meigen, 1830 mitochondrion equaling 6.2 × 10–8 mutations per site per generation (Haag-Liautard et al. 2008). The generation times of glacier-dwelling tardigrades, and of C. kaczmareki in particular, are unknown. However, to get rough insight, we used published data on another cryophilic tardigrade species Acutuncus antarcticus (Richters, 1904), where Tsujimoto et al. (2015) suggested 17 days generation time when cultivated at 15 °C. We stress however, that these numbers have to be taken with extreme caution as growth rates and generation times of ectothermic organisms scale with temperature [although not in a linear fashion, e.g., (Bartsch 2002)]. We thus presume that this estimate represents a minimum generation time for our focal taxon and that the real value for glacier species, like C. kaczmareki, would most likely be longer. Nevertheless, assuming an average growing season of 3 months, which roughly corresponds to the length of the season when the glacier surface is uncovered by snow, and the formation of cryoconite holes is possible both in high latitudes and altitudes, application of Tsujimoto et al.’s (2015) data would translate into 5.2 generations per year. We used this as an upper estimate to recalculate the divergence times from the tree made on 3rd positions. We note that such a calibration would downscale the time estimates, if the generation time is likely longer on glaciers. This is nonetheless desirable since we were primarily interested in assessing whether isolated C. kaczmareki populations represent a recent introduction via immigration after the last glacial maximum (LGM).

As a second approach, we calibrated the COI phylogenetic tree reconstructed from all positions with the estimates of mutation rates of the entire COI gene of northern polar arthropods calibrated by the divergence of the Beringian species, suggesting an average divergence rate of ~ 5.1% per million years (My) (Loeza-Quintana et al. 2019).

Phylogeographic analysis

The geographical component of the phylogeographic patterns was evaluated with the spatial analysis of molecular variance using SAMOVA version 1.0 (Dupanloup et al. 2002). It employs a simulated annealing procedure to cluster the sampling sites into a user-defined number of groups (K), taking the geographical locations into account, in order to maximize the proportion of the total genetic variance between groups (FCT) and minimize the proportion of variation among sites within groups (FSC). For this analysis, both sites in Himalaya (Baltoro and Chamser Kangri glaciers) with single successfully sequenced specimens were pooled together. We further excluded Nordenskiöldbreen from the analysis as it was represented by a single specimen.

To test whether gene flow exists between geographic regions, we jointly estimated the divergence times, migration connectivity, and sizes of analyzed populations using the Isolation with Migration MCMC algorithm (Hey 2010a, b). We originally set the analysis for three-population clusters assigned by SAMOVA: (1) Ebbabreen and Nordenskiöldbreen glaciers on Svalbard, (2) Qilian and Himalaya regions, and (3) Tien Shan Mts. As the COI gene phylogenetic analysis indicated that the Svalbard population is the most divergent, we performed a run using the [(1, 2), 3] population topology in IMa2 and did not integrate overall possible population trees, which is difficult with a single locus. After several short trial runs to assess priors and MCMC mixing, we performed two independent long runs, each with a 100,000 burn-in step and a 1,000,000 simulation step, sampling every 100th generation and swapping among 20 heated chains. Independent long and short trial runs showed almost identical results and were thus combined in a single analysis using the LogCombiner subprogram in the BEAST package. Since the entire COI fragment was used for this analysis, we used the aforementioned COI calibration by northern polar arthropods (Loeza-Quintana et al. 2019) to translate estimated parameters into absolute numbers.

As all runs indicated a very recent split between both Asian population clusters (see the Results section), suggesting little or very recent divergence among Asian samples, we also performed a two-population run assuming divergence between Svalbard (0) and Asia (1) only. To test for the possible long-term isolation between the analyzed populations as well as for differences in the long-term effective population sizes, we performed a Likelihood Ratio Test (LRT) of two alternative nested models using the L mode of IMa2. First, we kept all migration rates between Asia and Svalbard at zero, simulating vicariance and strict isolation between both regions. Second, we set both recent population sizes to equal values, thereby testing for unequal population sizes.

Demographic histories of C. kaczmareki and P. glacialis populations were estimated using Tajima’s D in DnaSP version 5.0 (Librado and Rozas 2009) and with the Extended Bayesian skyline plot analysis implemented in BEAST v.1.8.0 (Drummond et al. 2012), with subsequent analysis in Tracer v1.5 (Rambaut et al. 2009). In the case of C. kaczmareki, this analysis was restricted to individuals from Central Asia since all Svalbard samples possessed the single diverged haplotype and we, therefore, wanted to minimize the problem of population structure, which would tend to inflate coalescent times within the presumed population. The clock rate was set to strict, and the mutation model was set to Tamura-Nei with invariant sites site heterogeneity (as suggested by jModelTest) and base frequencies set to empirical values.

DNA metabarcoding of non-glacial habitats

In order to examine the possible presence of the putative glacial specialist C. kaczmareki in non-glacial habitats, we applied the metabarcoding approach to a set of environmental samples collected in the vicinity of Svalbard glaciers, including periglacial sediments and bryophyte samples from the foot of the glaciers and soil crust and tundra soil from the Petuniabukta area near Ebbabreen and Nordenskiöldbreen glaciers. A list of sampling sites is presented in Online Resource 3, and conditions for extraction of living animals, library preparation, and data processing leading to identification of denoised Zero-radius Operational Taxonomic Units (zOTU) are given in Online Resource 4.

Results

Cryoconicus kaczmareki and Cryoconicus antiarktos distribution

A review of the literature and inclusion of new material covered 60 glaciers worldwide (Fig. 2; Online Resource 1). The morphological examination and DNA analysis revealed several new localities for C. kaczmareki outside the reported range of Central Asia, i.e., Baltoro Glacier in the Karakoram, Kang Yatse, and Chamser Kangri glaciers in the Himalayas (Asia), and in the Central Spitsbergen (Svalbard Archipelago; High Arctic). Svalbard was examined in greater detail with resampling campaigns during 2014, 2016, 2017, and 2022, which provided several C. kaczmareki specimens. On Ebbabreen, we found 2 specimens from 6 cryoconite samples in 2014 identified as Ramazzottiidae in Zawierucha et al. (2016a); 13 specimens from 10 cryoconite samples in 2016; 14 specimens from 4 cryoconite holes in 2017; and 6 specimens from 4 cryoconite holes in 2022. A single specimen from the Nordenskiöldbreen was found in 2016 among a collection of ~ 5000 glacial tardigrades. In addition, the examination of DNA barcodes revealed the presence of a single specimen of C. antiarktos on the Ebbabreen in 2016, discovered only by DNA. Additionally, we found a population of C. antiarktos in the moss sample from Antarctica (Alexander Island: Abalation Valley) and used it in this study.

Morphological characters of the Svalbard C. kaczmareki specimens show differences from the Central Asian ones by the presence of clearly visible gibbosities in the hind legs (Fig. 3). Cryoconicus antiarktos from Ablation Valley differs from C. kaczmareki by presence of buccal tube armature (Fig. 4).

Fig. 3
figure 3

Cryoconicus kaczmareki from Svalbard: A Habitus, latero-ventral view; B Cuticular sculpture on the dorsal side; C Buccal apparatus, insert shows the enlargement of the tip region of mouth part to demonstrate the lack of teeth under PCM; D Fourth pair of legs, arrowheads indicate gibbosities; and E Leg without gibbosite. A, B Phase contrast microscope. C, E Differential interference contrast microscope. Scale bars in µm

Full size image
Fig. 4
figure 4

Cryococnicus antiarktos from Ablation Valley, Antarctica: A, B Habitus, latero-ventral view; C, D Claws. Arrowhead—claw base, asterix—cuticular bar at the claw base; E Buccal apparatus; and F Mouth part, arrowhead indicates oral cavity armature. Scale bars in µm

Full size image

Using previously published sequences and all the newly sequenced individuals of both Cryoconicus species, we aligned 459 bp of COI and 751-bp 28S rDNA. Newly obtained haplotypes were deposited into GenBank under accession numbers OP321127-OP321129 for 28S rDNA and OP316876-OP316879 for COI. Sequences of C. kaczmareki from Ebbabreen and Nordenskiöldbreen (Svalbard) were very similar to those from Central Asia sharing the same allele at the 28S rDNA locus, irrespective of their origin (Online Resource 5). At COI, we revealed 21 segregating sites (16 were parsimony informative), altogether defining 11 unique haplotypes across the complete C. kaczmareki dataset. The COI diversity was greatest within the Tien Shan population cluster, which appeared paraphyletic to individuals from both Qilian and Himalayan regions (Fig. 5a–c; Online Resource 5). One haplotype, central to the entire network, was present in all three Asian regions. The population from Svalbard (Ebbabreen and Nordenskiöldbreen) appeared the most divergent from other sites sharing a unique haplotype diverging by 6 mutations from the nearest Asian relative (1.307% p-distance).

Fig. 5
figure 5

a Ultrametric Bayesian phylogenetic trees of mtDNA COI variability of Cryoconicus kaczmareki reconstructed from neutral third codon position and b from all codon positions. 95% HPD intervals of node heights are provided by blue bars and bootstraps values are displayed above branches. Trees were rooted with C. antiarktos (not shown). Colored vertical bars indicate the samples taken in Svalbard (blue) and in Asia (red). c TCS haplotype network of C. kaczmareki COI dataset with colors corresponding to sampling sites and circle diameter to sampling sizes. d Plot of the relationship between genetic (ΦST) and geographic (km) distances between C. kaczmareki sampling sites. Red color indicates population pairs within Central Asia, blue denotes population pairs between Svalbard and Central Asian samples e TCS haplotype network of P. glacialis COI dataset with colors corresponding to sampling sites and circle diameter to sampling sizes

Full size image

In the Ebbabreen C. antiarktos sample, the 28S rDNA differed from the sequences of Antarctic individuals by 0.29% p-distance, being represented by four single-nucleotide polymorphisms (SNPs), two occurring in heterozygous states together with the shared variant and two at homozygous states (Online Resource 5). Three COI haplotypes were revealed altogether. The most dominant was shared by all C. antiarktos individuals from GenBank as well as by a newly sequenced individual from Ablation Valley, while two new singleton haplotypes differed by a single SNP (0.22% p-distance) from the dominant Antarctic haplotype. One belonged to the other Ablation Valley individual and the second to the Ebbabreen individual (Online Resource 5). Both new haplotypes were confirmed by reamplification and resequencing.

Intraspecific differentiation of C. kaczmareki and estimated divergence times

Phylogenetic analysis indicated that the Svalbard C. kaczmareki haplotype represents a sister lineage to Asian samples. The two molecular clock calibrations provided a rather broad range in estimated divergence times (Fig. 5a–b). Based on a neutral mtDNA mutation rate calculated from the 3rd positions and using the aforementioned lower estimate of generation time (5.2 generations per year), the median estimate for the split occurred ~ 99 thousand years ago (Kya) (95% CI: 56–155 Kya). On the other hand, when applying the COI gene calibration by Northern polar arthropods to the complete alignment, we found that the split occurred ~ 517 Kya (95% CI: 294–756 Kya).

The SAMOVA analysis indicates that geographical partitioning of C. kaczmareki mtDNA variability is best described by three- or four-group spatial structures that account for ~ 57% of variation attributable to among-group differentiation (Online Resource 6). The three-group scenario implies differentiation among Svalbard, Tien Shan, and populations located in Qilian together with the Himalayas, while the four-group scenario suggests an additional separation between Grigoriev Glacier and the remaining Tien Shan sites, thereby providing an uninformative grouping of single glaciers. Mantel test indicated a significant correlation between matrices of geographic and genetic ΦST distances among sampling sites (r = 0.941; p = 0.03), but this was merely due to the considerable differentiation of Asian samples from distant Svalbard population, while there was no evidence for isolation by distance (IBD) among Asian sites (Fig. 5d).

The Isolation with Migration (IM) analysis assuming three populations, as suggested by SAMOVA, i.e., (0) Svalbard, (1) Qilian & Himalayas, and (2) Tien Shan consistently indicated low migration rates among populations (probability mass peaks located near zero; Fig. 6, left column) and considerable differences in the effective population sizes with the largest population in Tien Shan and the smallest in Svalbard. Molecular clock calibration by Beringian divergence of Northern arthropods placed the maximum posterior probability distribution of the presumed split of Svalbard from Asia at ~ 310 Kya. The maximum posterior estimate for the split of both Asian population clusters was located at zero, suggesting their split either occurred very recently or we do not have enough signal in the data to detect the finer structure.

Fig. 6
figure 6

Result of Isolation with Migration (IM) analyses (right and left columns; a–f) and of Extended Bayesian Skyline Plot (g). Parameter values have been recalculated to absolute population sizes, split times in years, and connectivity in numbers of immigrants using the COXI rate for Arctic benthic invertebrates (see methods). Panels a, c, e: results of IM analysis assuming a three-population scenario, where population subscripts are as follows: 0—Svalbard; 1—Tien Shan; 2—Qilian & Himalaya; 3—the common ancestor of 1 and 2; and 4—common ancestor of all populations. a) Population sizes, c) split times, e) and migration rates. Panels b, d, and f: results of IM analysis assuming the two-population scenario, where population subscripts are as follows: 0—Svalbard; 1—Asia; and 2—the common ancestor of all populations. Results are reported for full isolation with migration (subscript “im”) as well as for pure isolation (subscript “isol”) models. b population sizes, d split times, f and migration rates. In order to merge different parameters into a single plot, the posteriors were rescaled with a scaling factor noted in the legend. Panel g Extended Bayesian skyline plot of cytochrome oxidase I marker for Cryoconicus kaczmareki Asian population (green) and Pilatobius glacialis middle Svalbard population (violet). The X-axis represents the time in mutation rate per site, and the Y-axis represents ln of Effective population size with solid lines depicting means and colors at the 95% highest posterior density (hpd) confidence intervals

Full size image

We also analyzed a two-population scenario assuming a single split between (0) Svalbard and (1) merged Asian samples (Fig. 6, right column). Negligible migration rates and considerable differences in population sizes between both regions were again suggested. In this analysis, the highest posterior estimate of their split time was at ~ 420 Kya, which is close to the aforementioned estimate from the sole phylogenetic tree, albeit the curve was relatively flat, suggesting that the data may not contain sufficient information to jointly estimate migration rates and split time. LRT of a nested model assuming zero migration rates did not reject the pure isolation model (2LLR = 4.037, df = 2, p value = 0.13), suggesting that the hypothesis of a gene flow between Asian and Svalbard populations has no significant support in our data. Both the full and strict isolation models otherwise provided similar estimates of population sizes and timing of splits (Fig. 6, right column). The LRT significantly rejected the model assuming equal population sizes between both regions (2LLR = 25.17, df = 1, p value < 0.001), corroborating the hypothesis that the Svalbard population has a lower effective size compared to Asian glaciers.

To get insight into the possible demographic fluctuations of C. kaczmareki, we measured Tajima’s D for the total dataset, Asian region only, and the Tien Shan population only, which possessed a reasonable sample size (Table 1). All values were close to zero and non-significant (total dataset D = 0.051; Asian dataset D = − 0.139, Tien Shan dataset D = 0.566). The reconstructed Extended Bayesian skyline plot suggested that the population size of the Asian region remained generally constant without notable changes during the reconstructed history (Fig. 6, middle panel).

Comparisons to dominating Arctic species Pilatobius glacialis

For comparative purposes, we evaluated the intraspecific diversity and demographic history of the Arctic glacier dominant P. glacialis, which was sampled from four glaciers of Central Svalbard, i.e., from the same biogeographical zone as the Arctic population of C. kaczmareki. COI haplotypes were deposited into GenBank under accession numbers OP316880-OP316885. All four samples of P. glacialis have higher molecular diversity than the local population of C. kaczmareki (Table 1 and Fig. 5c vs. 5e). Given the different sample sizes of both species, we further evaluated the probability of observing a single haplotype among exactly thirteen P. glacialis individuals, which corresponds to the sample size of C. kaczmareki on Ebbabreen. Since Arlequin software indicated no significant pairwise Φst differences between the four P. glacialis samples, we sampled with replacement 10,000 times exactly 13 individuals of P. glacialis from the total dataset and never observed a single haplotype among the simulated dataset. We also repeated the same procedure for each glacier separately, again finding a negligible probability of observing a single haplotype among 13 randomly selected samples (Table 1).

The P. glacialis dataset showed no evidence of major population fluctuations as Tajima’s D values from individual glaciers were non-significant (Table 1), and the EBSP indicate stable population sizes (Fig. 6, middle pane). A possible indication of some growth in the recent past was indicated by the mean of posterior values but not by the median.

Screening for Cryoconicus presence outside glaciers from a metabarcoding dataset

Metabarcoding screening of periglacial habitats delivered 438,481 R1 reads after the bioinformatic filtering steps (average number of reads per sample: 33,729). Raw sequencing reads are openly available in SRA at https://www.ncbi.nlm.nih.gov/sra, with accession numbers available in Online Resource 3. The analysis revealed 807 zOTUs, of which 63 had 95% or higher sequence similarity hit to published tardigrade sequences, including P. recamieri (a sister species of P. glacialis) reported to dwell in polar tundra (Zawierucha et al. 2016b; Gąsiorek et al. 2017) that was detected in abundant read numbers in a coastal soil sample from Southern Svalbard. However, no zOTU matched the published or newly acquired sequence of C. kaczmareki or C. antiarktos, suggesting these species were not recovered in our dataset. Similarly, BLASTn search of Cryoconicus and P. glacialis haplotypes against all quality-filtered reads revealed no match to any of these taxa in any sample.

The positive control of mixed gDNA revealed 70,755 retained reads, of which 25% were assigned with > 99% identity to C. kaczmareki Svalbard haplotype and 29% to P. glacialis. This corroborated our approach to detect the reads of these glacial specialists among the mixture of several other tardigrade and rotifer species.

Discussion

Apparent vs. real discontinuities in a global scale distribution of Cryoconicus

In this study, C. kaczmareki and C. antiarktos have been recorded in small numbers at two High Arctic glaciers (Ebbabreen and Nordenskiöldbreen), far outside their terra typica. Arguably, such an observation could not be an artifact of PCR contamination since the analyses from Victoria Land, Antarctica (Guidetti et al. 2019), Central Asia (present study and Zawierucha et al. 2018b), and Maritime Antarctica (present study) samples were conducted in different laboratories or different years (and different chemicals), and the COI and 28S haplotypes from the Arctic (Svalbard) and Antarctic (Alexander Island: Ablation Valley) were resequenced by an independent sequencing service to eliminate potential errors. Differences in haplotypes between Arctic and Asian samples of C. kaczmareki (Fig. 5) further argue against human-mediated contamination between Asia and Svalbard. Both species thus appear to have transcontinental distribution but with apparent large-scale discontinuities.

The most trivial explanation of such disjunct populations would be the unreported presence of both species in intermediate areas. However, large body size and densely pigmented (blackish/brownish) color makes the Cryoconicus genus easily distinguishable from other tardigrades and, therefore, unlikely to have been overlooked in previous studies, which investigated literally thousands of individuals collected through various seasons from many worldwide glaciers (Fig. 2 and the Online Resource 1). Other black tardigrades inhabit glaciers, for instance, the genus Cryobiotus in the Alps and the Himalayas or the genus Kopakaius in New Zealand (Zawierucha et al. 2023), but they differ in shape and claw structure, which are visible even under a stereomicroscope. Hence, although we cannot exclude the possible occurrence of these Cryoconicus species in some uninvestigated areas (Northern Andes, West-Pacific Coast of USA, Northern Greenland, Ural, the Caucasus), it appears their distributions extend across several major biogeographical zones and contain large gaps spanning thousands of km.

The research on the biogeography of small organisms has followed the ‘everything is everywhere’ theory, postulating that species have large or even cosmopolitan ranges, but a recent merging of traditional and DNA taxonomy has changed this view, showing that even small organisms have clear biogeographical patterns with rarely global distributions (Green et al. 2004; Fontaneto et al. 2008; Wu et al. 2011; Porco et al. 2012; Bates et al. 2013; Iakovenko et al. 2015; Fontaneto 2019; Morek et al. 2021). Cryoconicus antiarktos thus appears as one of few animal species with truly bipolar distribution (see Kaczmarek et al. (2020) for Paramacrobiotus fairbanksi Schill et al., 2010) and adds to the limited number of tardigrade species diagnosed by combined DNA and morphological analyses that are distributed in more than one zoogeographic zone (Stec et al. 2020; Morek et al. 2021; Guil et al. 2022).

Only a single specimen of C. antiarktos was extracted from a sample from Ebbabreen and it differed by a single SNP from mtDNA haplotypes reported from the Antarctic populations. Hence, we cannot infer whether there is a viable population in the Arctic. In comparison, C. kaczmareki has been repeatedly encountered on the Ebbabreen during repeated sampling over four years. The population genetic analyses did not indicate any recent gene flow between the Asian and Arctic populations, suggesting they have been isolated in both historical as well as contemporary periods. The split time between both populations was estimated within a large range, though. The estimate for entire COI gene using the molecular clock of Arctic benthic invertebrates indicated an ancient split at 294–756 Kya (or ~ 420 Kya when jointly estimating isolation & migration parameters with IMa), while applying neutral mtDNA mutation rate to 3rd COI positions with 5.2 generations per year indicated more recent split at 56–155 Kya. It is likely that mutation rates and generation turnover in cold environments are lower than in temperate conditions and so it is possible that the latter estimate might have been biased toward lower values, while the real split is older. Indeed, the lack of relevant calibration points makes the application of molecular clock challenging for our taxon. Nonetheless, both methods, albeit relying on different rationale and approaches, agreed on the assumption that the split time between Asian and Arctic populations predates LGM, thereby supporting the hypothesis that C. kaczmareki established a long-term stable and genetically distinct population in the Arctic.

Mechanisms causing the disjunctive ranges

The next explanation for such disjoint ranges of glacial tardigrades would be that the northern glacial sheets were much greater during ice ages (Sejrup et al. 2005; Hughes and Gibbard 2018) and currently isolated Arctic population would thus represent fragmented relics of once widespread populations interconnected to other regions by short-range migration among now-extinct demes. However, the ice sheet connection between Central Asia and the Svalbard Archipelago has not existed even during the greatest Pleistocene glaciations. Thus, even if C. kaczmareki might have occurred on potential stepping stones (such as expanded glaciers in Turkey, the Urals or Novaya Zemlya Archipelago), its range would still include large gaps spanning hundreds or thousands of kilometers. This, in turn, indicates long-distance dispersal that generated the disjunct ranges of both our species.

While microinvertebrates are traditionally supposed as good dispersers (Bertolani et al. 1990; Mogle et al. 2018; Fontaneto 2019), the mechanisms of how tardigrades disperse are poorly known (Jørgensen et al. 2007; Gąsiorek et al. 2018; Morek et al. 2021; Guil et al. 2022). One dispersion vector commonly assumed for microinvertebrates is the wind (Faurby et al. 2008; Dabert et al. 2015; Ptatscheck et al. 2018; Zawierucha et al. 2019a) but relatively little is known about the importance of anemochory for meiofaunal dispersal across large scales (Jørgensen et al. 2007; Dabert et al. 2015; Fontaneto 2019; Rivas et al. 2019). The prevailing Holocene Arctic wind flow has been a westerly circumpolar flow meeting northeasterly trade winds on the southern edge (e.g., Barry and Chorley 2009), and recent data suggest that Arctic dust deposits may have been transported from Asia during the LGM (Újvári et al. 2022). This may theoretically account for C. kaczmareki distribution. However, it is unlikely that the isolating circumpolar winds of Antarctica and several contrary wind patterns between the North and South poles would permit bipolar distribution of C. antiarktos via anemochory.

Another mean of dispersal may involve migrating birds, which have been suggested to disperse North American glacier ice worms (Dial et al. 2016a; Hotaling et al. 2019) and to shape the spatial patterns of central American tardigrades during the Great American Biotic Interchange (Mogle et al. 2018). Since migratory birds may rest on glaciers and snow (Rosvold 2016), it is tempting to speculate about the transport of C. antiarktos to Svalbard, e.g., by Sterna paradisaea Pontoppidan, 1763, which is famous for its yearly bipolar migration (Egevang et al. 2010). On the other hand, while some bird species do migrate between Central Asia and the High Arctic (Snell et al. 2018), the C. kaczmareki arctic population diverged from Asia long ago, and the historical bird migratory routes are unknown and might have substantially changed during Pleistocene glaciations (Zink and Gardner 2017; Somveille et al. 2020).

Implications for population stability of (extremely) rare taxa

The comparison of the isolated C. kaczmareki Arctic population with other Arctic glacier-dwelling tardigrade species and with the population on Asian glaciers has interesting connotations for two important questions of general ecology: (i) whether the species abundances remain stable in time and (ii) whether species are suffusively or diffusively rare through time, i.e., rare throughout a time period or rare at some times, but not at others (Colinvaux 1979; Gaston 1994).

Specifically, while C. kaczmareki is a common species on Asian glaciers with a large effective population size inferred from mtDNA variability (Fig. 5c; Zawierucha et al. (2018b)), the isolated Arctic populations persist in relatively low abundances, being detected only on Ebbabreen and Nordenskiöldbreen. The species has not been found anywhere else across dozens of glaciers and surrounding habitats in the Svalbard archipelago. Indeed, it is possible that this species may persist in higher numbers in other unsampled habitats, e.g., the accumulation zones of glaciers, which may potentially host large populations of tardigrades (Shain et al. 2021) and our samples may represent only the population sink from such nearby habitats. However, the low genetic variability of the Ebbabreen and Nordenskiöldbreen population (single COI haplotype) suggests that C. kaczmareki has a significantly lower effective population size in the Arctic than in Asia, making it unlikely that sampling sites on Ebbabreen and Nordenskiöldbreen are replenished from any nearby (unsampled) reservoirs.

Cryoconicus kaczmareki is thus much rarer in comparison to other glacier Eutardigrada (P. glacialis) that are encountered in densities orders of magnitude higher in the same sites (Zawierucha et al. 2016a, 2019b, 2020; Novotná Jaroměřská et al. 2021). Such differences in relative abundances have likely persisted throughout long historical periods since currently dominating P. glacialis has a significantly higher genetic diversity on glaciers of central Svalbard (Table 1), suggesting its larger long-term effective population size compared to C. kaczmareki. Moreover, we found no evidence of demographic fluctuations in C. kaczmareki and P. glacialis core ranges (Asia and Central Svalbard, respectively), suggesting that the current dominance of P. glacialis on Svalbard glaciers is not a result of a recent population expansion. Hence, while C. kaczmareki forms a dominant taxon on Asian glaciers, it seems to constitute a suffusively rare component of the Arctic glacial ecosystems.

While rare species are threatened by ecological drift and stochastic processes like local fluctuations to a much greater extent than more common species (McKinney et al. 1996; Kunin and Gaston 1997), they may mitigate increased extinction risk by several mechanisms (rev. in Kunin and Gaston 1997) such as adaptation to a broader niche (eurytopy) (e.g., Futuyma and Moreno 1988), occurrence in productive environments (Vermeij 1986), or wider distribution that buffers local extinctions (Maurer and Nott 2001). However, none of these premises are imminently applicable to the long-term endurance of the isolated population of C. kaczmareki in the Arctic. The isolated Arctic populations do not seem to be part of a larger metapopulation, making it unlikely to mitigate the increased extinction risk of local demes. The possibility of eurytopy is also put into question by the fact that glacial habitats are generally considered nutrient poor (Porazinska et al. 2004; Zawierucha et al. 2018b; Poniecka et al. 2020), albeit the content of organic matter may sometimes be relatively high (Rozwalak et al. 2022).

Conclusion

Although microscopic animals play key roles in ecosystems, they are not as well studied as single-cell organisms. Our large-scale study corroborates growing evidence (e.g., Morek et al. 2021) that microscopic animals do not have truly cosmopolitan ranges. However, our findings, together with some earlier reports (Kaczmarek et al. 2020; Stec et al. 2020; Morek et al. 2021; Guil et al. 2022), suggest that interpolar- or at least intercontinental-scale movements of meiofauna, including cryophiles, do occur, especially when such organisms can rapidly establish new populations via asexual reproduction (Kokko and López-Sepulcre 2006; Fontaneto 2019). Cryoconicus kaczmareki distribution patterns further demonstrate that microinvertebrates may persist long-term in relatively low abundance and in a geographically extremely limited area. This is true even in unstable habitats, like on an Arctic glacier surface, where population density may extremely fluctuate from one day to another (Zawierucha et al. 2019b). Similar examples of locally persistent distributions have recently been reported in European Alpine plants (Smyčka 2019), suggesting that the inhabitants of extreme environments might often persist long term in isolated and extremely fragmented populations.

An increased focus on long-term and fine-scale monitoring of local populations of microinvertebrates may thus considerably impact our understanding of how ecosystems are assembled, which is mostly derived from macroscopic organisms. As a relatively easily distinguishable habitat with stable seasonal thermal regimes and truncated food webs, glaciers may prove to be excellent models and natural laboratories in this quest.