Intra‐Alpine Islands: Population genomic inference reveals high degree of isolation between freshwater spring habitats

Alpine spring ecosystems have long been considered as highly isolated, island‐like habitats. This presumption, however, is insufficiently supported empirically and conclusions about spring isolation have been based on indirect evidence. Therefore, we investigated the population genomic structure of Partnunia steinmanni Walter, 1906, a strictly spring‐dwelling water mite (Hydrachnidia) species, to shed light on the degree of interconnection among freshwater spring habitats.

In the sense of classic island biogeography theory (MacArthur & Wilson, 1967), increased habitat isolation is predicted to negatively impact immigration, and accordingly the re-colonization rate, thus fostering genetic bottlenecks (Broquet et al., 2010;Pinheiro et al., 2017). Consequently, reduction of populations' resilience through diminished abilities to adapt to changing environmental conditions becomes likely (Sgrò et al., 2011). Decreasing population size due to anthropogenic impact may further exacerbate this effect (Elsen & Tingley, 2015;Shama et al., 2011).
However, conclusions about the degree of isolation of spring habitats have mainly been based on assumptions derived from community composition changes within small-scale study areas and thus are primarily indirect evidence for the interconnection between habitat patches (e.g. Fattorini et al., 2016;Von Fumetti & Blattner, 2017). Studies explicitly investigating spring population interconnection based on genetic structure are still scarce, limited in marker resolution to few traditional genetic loci, and/or conducted in extreme environments such as deserts dominated by endorheic and subsurface basins that are difficult to compare to other environments (e.g. Adams et al., 2018;Myers et al., 2001;Stutz et al., 2010).
To empirically assess alpine crenic habitat interconnection, we investigated the crenobiontic water mite species Partnunia steinmanni Walter, 1906 (Figure 1) that exhibits a primarily alpine distribution area with additional records from the Tatra and Western European lower mountain ranges (Gerecke, 1993). The genus Partnunia Piersig, 1896 includes in total ten species described from Europe and Asia (Gerecke, 1996) and belongs to the Hydryphantoidea, a phylogenetically basal water mite superfamily that was recently recognized as monophylum (Blattner et al., 2019;Dabert et al., 2016;Di Sabatino et al., 2010). In addition to P. steinmanni, only the species P. angusta (Koenike, 1893) has been described in Central Europe and assigned to the genus Partnunia, which exhibits relatively homogenous morphology (Di Sabatino et al., 2010;Gerecke, 1993). Partnunia angusta is restricted to the Alps and northern Prealps and in addition to springs also appears in spring brooks and low order streams (Di Sabatino et al., 2010). We focused on the strictly crenobiontic P.
steinmanni that shows strong habitat preference for shaded springs dominated by moss and gravel substrate (Gerecke, 1993). Partnunia steinmanni were sampled in 12 springs located in different major protected areas across the Alps, and restriction site-associated DNA F I G U R E 1 Living Partnunia steinmanni specimen 1 mm sequencing (RADseq) was performed. Subsequently, we inferred intra-and inter-population genomic structure, determined the influence of spatial structure on the genetic differentiation, and calculated demographic estimates to assess the degree of spring habitat isolation and estimate the genetic diversity of a characteristic crenobiontic species.
2 | ME THODS 2.1 | Study sites, sampling and pre-processing Partnunia steinmanni specimens ( Figure 1) were sampled with a hand net (100 µm) in springs located in six main protected areas across the Alps in summer 2019 and 2020 (Table 1) (Table 1). Sampling permissions were granted for each protected area by the respective authorities.
Each spring population consisted of ≥30 individuals sorted out alive directly in the field, and each individual was transferred to a single well of a 48-well cell culture plate (Sarstedt AG & Co. KG, Nümbrecht, Germany) prefilled with water directly from the spring.
The plates were held at 4°C for one week until the specimens were subsequently transferred to molecular grade ethanol (100%) and stored at −20°C until further processing. This procedure resulted in a starvation period with the aim to reduce possible sample contamination due to residual gut content, which potentially can be detected until one week after starvation (Martin et al., 2015).

| DNA extraction, RADseq library preparation and sequencing
Each mite was first submerged in molecular grade water to remove tracted DNA, were chosen to be further processed, resulting in a total of 256 P. steinmanni specimens (Appendix S1).
Due to relatively low initial DNA yield (mean ± SD: 0.2 ± 0.15 ng/ µl, Appendix S1) that was obtained from the water mites, a multiple displacement amplification (MDA) was performed for each individual to increase the amount of DNA. This procedure has proven to work well for RADseq, without introducing notable genotyping bias as shown in previous studies (Blair et al., 2015;Cruaud et al., 2018;De Medeiros & Farrell, 2018). For MDA, the REPLI-g Mini kit (Qiagen, Hilden, Germany) was applied to each sample with 5 µl of DNA eluate resulting in 161 ± 23.7 ng (mean ± SD, Appendix S1) amplified genomic DNA.
To obtain individual-specific datasets, demultiplexing and quality filtering of the raw data was done with process_radtags, a component of the Stacks 2 V2.59 pipeline (Rochette et al., 2019), allowing one mismatch position in the barcode sequence (Appendix S1) and requiring an intact PstI restriction residual. Due to the use of bluntend ligated sequencing adaptors, resulting in mixed orientation of forward and reverse reads, the bestrad option was enabled in pro-cess_radtags that automatically scans for the barcode sequence on either read and corrects the orientation.
Following Lucek et al., 2020, ustacks V2.59 (Rochette et al., 2019 was used to de novo assemble the pre-processed restriction site flanking forward reads with a minimum stack size of 50 reads, three allowed mismatches between reads to be associated with the same stack, disabled gapped-alignment and disabled haplotype calling from secondary reads. Reverse reads were not included because of high read overlap due to relatively small final insert sizes (± 200 bp) and therefore sequence redundancy. To account for putatively high interpopulation divergence caused by the large geographic extent of the study, the assembly was performed separately for individuals belonging to a specific sampling area resulting in six different assemblies. Homologous RAD loci were then identified between the different assemblies using swarm V3.1 (Mahé et al., 2021), merging contigs with 97% sequence similarity. Only contigs that occurred in all six initial assemblies were retained.
Resulting contigs were subsequently filtered for contaminants by blasting each contig against the NCBI GenBank nucleotide database (Agarwala et al., 2018(Agarwala et al., , accessed: 15.3.2021 with the blastn function as implemented in BLAST+ V2.2.23 (Camacho et al., 2009). This yielded 1516 final contigs that were not identified as derived from extraneous origin (e.g. bacteria, archaea or viruses).
The pre-processed raw reads were then aligned against the de novo assembly using minimap2 V2.17 (Li, 2018) in short read mode, and the resulting alignments were converted, sorted, and indexed with SAMtools V1.12 (Danecek et al., 2021). Next, the mpileup und call functions of BCFtools V1.10.2 (Danecek et al., 2021) were used with a maximum read depth filter of 5000 to call a total of 7069 SNPs. This preliminary genotype dataset was then filtered with VCFtools V0.1.17 (Danecek et al., 2011). Sites with more than 60% missing data, minor allele frequency of <0.01 and minor allele count of <3, genotype quality of <Q20, and read depth of <5× and >800× at the individual level were removed, as well as indel positions. Furthermore, 18 individuals showing >56% missing data were removed. This approach resulted in a final dataset of 2263 SNPs for 238 individuals with a mean read depth of 172× that was used for all subsequent genomic analysis.

| Population genomic and phylogenetic inferences
Physical variant linkage pruning was conducted in PLINK V1.90b6.24 (Chang et al., 2015) with a sliding window size of 50 variants, a step size to shift the window of 10 variants, and a r 2 threshold of .1 to obtain a SNP dataset with a subset of markers that are in approximate linkage equilibrium. Linkage pruning resulted in 1409 SNPs that were used to infer population structure with Admixture V1.3 (Alexander et al., 2009). Admixture was run assuming between 1 and 20 putative populations (K) and the most probable K was assessed by the Admixture cross-validation procedure. Additionally, a principal component analysis (PCA) was calculated in PLINK for the linkage pruned markers. Both analyses were post-processed and plotted in R V4.1.0 (R Core Team, 2021).
The relationship among all individuals was assessed with RAxML V8.2.12 (Stamatakis, 2014) under a generalized time-reversible (GTR) model of evolution with optimized substitution rates and gamma model of rate heterogeneity on the dataset without linkage pruning. To avoid disconcerting influence of admixed individuals on the general tree topology (Seehausen, 2004;Shirk et al., 2021), individuals showing >20% admixture were excluded. Because only polymorphic sites were used, a Lewis ascertainment bias correction (ASC_GTRGAMMA function in RAxML) was implemented, and statistical branch support was evaluated by computing 5000 bootstrap replicates. The resulting unrooted bipartition tree was visualized in TreeViewer V1.2.2 (Bianchini, 2021).
Isolation by distance (IBD) was assessed by first computing geographic distances in km between the sampling location as implemented in the geodist V0.0.7 (Padgham & Sumner, 2021) R package (Appendix S2). GenoDive V3.05 (Meirmans, 2020) was then used to calculate pairwise F ST values between all combinations of P. steinmanni spring populations (Appendix S3), followed by a Mantel test (Mantel, 1967) with 20,000 permutations between log transformed geographic and genetic distance matrices. To further quantify the amount of genetic variance explained by the spatial structure, we performed a distance-based redundancy analysis (dbRDA) as implemented in the vegan V2.5-7 R package (Oksanen et al., 2020) (Table 1) and MEMs were then included in the dbRDA as explanatory variables to represent spatial structure.
Furthermore, an analysis of molecular variance (AMOVA) (Excoffier et al., 1992) was calculated with 10,000 permutations based on pairwise F ST in GenoDive to assess the level of population differentiation. We tested for significant genetic differentiation among springs within, as well as between different sampling areas.

| Population genomic structure
The best-supported number of genetic clusters (i.e. lowest standard error of cross-validation error estimate) identified by Admixture was 12 (Appendix S4), corresponding precisely to the number of sampled springs (Figures 2 and 3a). Overall, each spring thus consisted of individuals exhibiting a spring-specific genotype. Furthermore, admixture tended to be slightly more pronounced between springs within sampling areas, for example between the VA4 and VF3 populations located in the Swiss NP area or between GSC and KOE in the Gesäuse NP (Figure 2).  (Figures 2 and 3c). Furthermore, the sprigs located in the central Alps, Swiss-and Adamello-Brenta NP, seem to be closer related to each other than to the peripheral sites, although statistical bootstrap support was weak.

| Population differentiation and isolation by distance
The analysis of molecular variance (AMOVA) revealed the strongest genetic differentiation among the sampling areas, that is, among national parks (29% of observed variation, F ST ± SD = 0.34 ± 0.005, p < .001). A significant fraction of the total variation, however, was also explained by the springs within the areas (21% of observed variation, F ST ± SD = 0.21 ± 0.008, p < .001). Pairwise F ST between the geographically closest spring populations within the Gesäuse NP, Swiss NP and the Adamello-Brenta NP was consistently higher than 0.18 (Appendix S3). In springs by spatial structure, that is, their geographic proximity and catchment areas (Figure 4).

| Genetic diversity and demographic estimates
Overall, relatively low and homogenous levels of genetic diversity was observed between springs (mean H e ± SD = 0.113 ± 0.023; mean H o ± SD = 0.122 ± 0.0273) ( Table 2). The RIF spring population showed the lowest expected and observed heterozygosity.
Interestingly, RIF and KOB also exhibit increased levels of inbreeding compared to the other populations. BRE, VF3, HIB and GSC showed no evidence for putative inbreeding. Effective population size for P. steinmanni was low and on average estimated at 9.9 ± 4.3 (mean ± SD) individuals (Table 2 and Appendix S6). Due to infinite confidence intervals, KOB failed at providing reliable N e estimates, likely caused by sampling bias due to the low amount of individuals processed (Do et al., 2014). The assessment of contemporary migration revealed very low migration rates (m) between (mean m ± SD = 0.01 ± 0.006) and genetic exchange almost exclusively

F I G U R E 4
Distance-based redundancy analysis (dbRDA) of genetic distances and spatial structure. Pairwise F ST were used as dependent variables and geographic distances (Appendix S2) as well as the main catchment areas across the Alps were included as explanatory variables to assess the influence of spatial structure on the P. steinmanni population differentiation. The first two plotted axis explain 31.  within springs (mean m ± SD = 0.87 ± 0.02). Only from VF3 to VA4, both springs located in the Swiss NP, a slightly higher migration rate of approximately 6% (m = 0.06; 1.2 out of 21 individuals) occurred ( Figure 5 and Appendix S7).

| DISCUSS ION
Spring ecosystems have been considered as putatively isolated insular habitats despite limited empirical evidence (e.g. Cantonati et al., 2006;Fattorini et al., 2016). To evaluate this concept, we here investigated the population genetic structure of a strictly spring-dwelling water mite species within and among protected areas across the Alps.
We show that even geographically close populations of P. steinmanni exhibit spring-specific genotypes, and vast genetic differentiation combined with limited inter-spring migration. This provides strong evidence of restricted dispersal and gene flow, consistent with an island-like habitat character of alpine spring ecosystems. Low degree of habitat interconnection between springs has already been inferred indirectly based on macroinvertebrate and spring-related stygofauna community composition data (Fattorini et al., 2016;Von Fumetti & Blattner, 2017), springs from extreme environments (Myers et al., 2001), or species that are not exclusively spring-bound (Engelhardt et al., 2011). Our study now offers population genomic support of the idea that alpine spring ecosystems represent islandlike habitats.
Landscape-dependent population structure is, among other factors, influenced by species-specific dispersal abilities that determine the relevance of putative barriers, such as topographic and environmental gradients (Garant et al., 2007;Storfer et al., 2010; Van Buskirk & Jansen van Rensburg, 2020). Springs are known to harbour diverse species assemblages with different dispersal capacities (e.g. Stevens et al., 2021), including strongly or weakly dispersing taxa (De Bie et al., 2012). Consequently, depending on the study organism, the degree of geographic isolation may differ, potentially leading to opposing conclusions about spring interconnection.
In contrast to other spring-dwelling taxa, water mites such as P.
steinmanni show live-stage-dependent dispersal abilities. As larva, they parasitize insect imagines with differing flight abilities, putatively allowing for effective dispersal, and are as adults restricted to a single spring (Martin & Stur, 2006;Zawal, 2003). The host species is, however, potentially impacted by parasite load, which may alter its flight capacities and consequently restricts its dispersal abilities and migration distances (Sánchez et al., 2015;Smith, 1988). Water mites F I G U R E 5 Contemporary migration rate estimates between spring populations. Directional migration was calculated by using BA3-SNP and is represented as arrows, whose sizes are relative to the proportion of migrating individuals and colour-coded according to their population of origin. Exact migration rates, and accordingly number of individuals relative to absolute population sizes can be seen in Appendix S7  (Engelhardt et al., 2011) and were shown across the Alps for both, terrestrial and aquatic taxa, to be associated with post-glacial recolonization from distinct refugia (Asztalos et al., 2021;Hewitt, 2000;Lucek et al., 2020). In contrast to that, the Swiss National Park exhibits a unique genomic setting with only a single individual in the VS6 spring showing some genetic similarity with the KOE spring of the Gesäuse NP region.
This may suggest that local topography could have acted as a barrier and hindered the western lineage putatively originating in the Rhône drainage basin from spreading directly northwards; however, denser sampling of possible contact zones would be necessary to support this assumption.
Apart from these general patterns, the overall high degree of genetic differentiation between P. steinmanni populations from different springs (Appendix S3-F ST ) may have originated as a result of many different intra-alpine glacial refugia possibly associated with main catchment areas, where the species was able to survive locally during glacial periods and experienced post-glacial isolation due to limited dispersal capacities, analogous to the alpine caddisfly species Drusus discolor (Rambur, 1842) (Pauls et al., 2006).
The observed genetic differentiation among springs could also reflect isolation caused by local adaptation to different microhabitats in springs. However, the influence of isolation by distance vs.
isolation by environment on population structure is discussed controversially (e.g. Aguillon et al., 2017;Sexton et al., 2014;Van Buskirk & Jansen van Rensburg, 2020) and seems to be taxon, environment and time-scale-dependent. By assessing the correlation between geographic and genetic distance, we were able to show the presence of isolation by distance (IBD) among springs. Partnunia steinmanni is known to be strictly restricted to alpine crenic habitats and prioritizes a specific benthic microhabitat, that is, gravel-dominated substrate rich in moss in forest or shaded locations (Gerecke et al., 2005(Gerecke et al., , 2009 (Martin et al., 2009). Isolation by population-specific host species preference could thus have further shaped the population structure of P. steinmanni. High degree of association with spring-specific host species assemblages could potentially lead to a restricted dispersal to springs at close proximity and consequently foster genetic isolation. The exceptionally low between-spring migration rates shown by our results strengthen this assumption. Indeed, a slightly higher migration rate occurred between VF3 and VA4, where there is a lack of topographic barriers such as mountain massifs between these two springs ( Figure 2).
However, to describe additional processes potentially causing the strong population distinctiveness and isolation of springs and assess the influence of IBE, further investigations and thorough evaluation of putative environmental differences between P. steinmanni habitats need to be conducted.
Effective population sizes (N e ) were generally low, that is around 10 individuals per population, which is comparable to other parasitic Acari (see e.g. Huber et al., 2019). Compared to the census population sizes that can easily exceed tens and even hundreds of P.
steinmanni individuals (e.g. Gerecke et al., 2009;Kreiner et al., 2018), the estimated N e is relatively low. This may indicate a reduction of effective population size due to low migration and high isolation over time, implying a bottleneck scenario with putatively increasing influence of genetic drift. Consequently, genetic diversity loss and the possibility of reduced fitness of the species may occur (see e.g. Broquet et al., 2010;Charlesworth, 2009;Hohenlohe et al., 2021).
Due to re-colonization of the periglacial inhabitable areas, we assume that founder effects could also have induced the relatively low N e and reduced genetic diversity within springs (see e.g. Montero-Pau et al., 2018;Peter & Slatkin, 2015). The effect of strong genetic differentiation among and low genetic diversity within populations particularly applies for species inhabiting formerly glaciated areas (Galbreath & Cook, 2004;Pečnerová et al., 2017) and has been shown for other water mite species (Bohonak, 1999). Suitable estimation of effective population size, however, can be influenced by sampling strategy (Do et al., 2014;Hare et al., 2011). Especially in species that show pronounced spatial structure and isolation, single sample strategies, as N e estimation based on LD, should be validated by temporal methods including samples of multiple generations (Neel et al., 2013).
The RIF spring population near Zermatt showed exceptionally low genetic diversity, and a relatively high inbreeding coefficient compared to the other springs, that revealed rather homogenous heterozygosity estimates. We assume that the geographic location of Zermatt, surrounded by pronounced mountain massifs, and therefore, high degree of topographic isolation (see Figure 2) has potentially driven this pattern. The comparably high N e that was estimated for this population may result from the above-mentioned sampling bias and should be validated by increasing the sample size, that is number of individuals.
Taxonomically, the description of P. steinmanni has been based on few individuals from very different geographic locations, including the locus typicus in the Rifelbord Wildlife Area near Zermatt, but also far distant springs in eastern Austria and even outside the Alps (Gerecke, 1993). The morphological characters of the a priori morpho-species are considered to be rather variable; thus, the ap- propriateness of the species delimitation should be questioned (R. Gerecke, personal communication). Furthermore, intra-species clade separation has been shown for P. steinmanni by investigating traditional genetic species delimitation markers (Blattner et al., 2019).
Given the strong population structure, a thorough re-evaluation of the P. steinmanni morphospecies should be conducted to assess the possibility of P. steinmanni being a species complex rather than a single, well-defined species.
To conclude, our results provide strong evidence for a high degree of insularity of alpine spring habitats, likely shaped by Pleistocene isolation in different intra-alpine refugia associated with main catchment areas as previously shown to be characteristic for alpine headwater environments. Crenobiontic species with limited dispersal abilities and low inter-population migration rates such as P. steinmanni show restricted gene flow, resulting in high genetic diversity and the potential for speciation across the Alps. Each spring and study area contains populations with a unique genetic make-up, and the loss of individual habitats will directly drive a decrease in overall genetic diversity. Our work thus highlights the importance of protected areas, such as National Parks that limit anthropogenic impact, as archives of genetic biodiversity.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interests to declare.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13461.

DATA AVA I L A B I L I T Y S TAT E M E N T
Short read data were uploaded to the European Nucleotide Archive (ENA) and are available under the study accession number: PRJEB47010.

B I OS K E TCH
Lucas Blattner is working in the field of molecular biogeography and is interested in environmental sciences in general. This study is part of his Ph.D. thesis, which focused on studying alpine spring ecosystems by investigating crenobiontic Hydrachnidia species to understand their ecology and mechanisms shaping their distribution.
Author contributions: All authors edited and approved the manuscript. LB conceived and designed the study and wrote the manuscript. He planned and conducted field and laboratory work as well as bioinformatic processing of the data. KL substantially contributed to data analysis, interpretation and manuscript editing.
NB was mainly involved in laboratory protocol development and optimization. DB supported the study with conceptual and data analysis input and SF contributed to conceptualization, manuscript editing and study realization.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found in the online version of the article at the publisher's website.