East Meets West in the Western Carpathians: Contact Zone Between Alpine and Carpathian Caddisy Lineages Evidenced by COI Barcodes

European mountains are important areas regarding biodiversity of the continent, and they also harbour diverse freshwater fauna, which is critically endangered in terms of the current species extinction. However, sucient knowledge of this valuable part of European biota is no longer possible without molecular data. This study focuses on the genetic diversity and distribution patterns of the classical representative of the mountain freshwater fauna, caddisy Rhyacophila tristis, in the Western Carpathians. Based on the COI mitochondrial marker, two genetic lineages (separate BINs) were identied. BIN_W was found in 16 localities in the western part of the study area, BIN_E in 44 eastern localities. The data obtained indicate that BIN_W occurs in a signicantly narrower altitude range, BIN_E is more closely related to the Balkan populations than to co-occurring BIN_W, and that the contact zone of the lineages passes through the Western Carpathians. The study revealed phylogeographic and demographic differences between lineages, supporting hypothesis of their evolutionary independence and specic ecological preferences. The obtained genetic data shifted our knowledge on the populations of the studied caddisy and suggested patterns that could be common to other freshwater species. This could help us to protect unique freshwater ecosystems and preserve European biodiversity.

The knowledge on the spatial distribution of genetic diversity at a species level is an essential prerequisite for understanding the resilience of biota in the changing environment across geological periods (Lande & Shannon, 1996;Frankham, Briscoe & Ballou, 2002). This is especially important for species inhabiting speci c biotopes such as mountain aquatic ecosystems. The caddis y Rhyacophila tristis Pictet, 1834 is one of the three species from the Rhyacophila tristis group, widespread in Europe. Adults of R. tristis are relatively good iers, which guarantees e cient dispersal. They prefer upper parts of streams, because larvae require well-oxygenated water with low levels of organic matter (Robert & Curtean-Bănăduc, 2005). Bálint et al. (2011) provided the rst data on the large-scale population genetics of this (sub)montane species on the European level. The study revealed strong genetic differences between populations from the western and eastern part of Europe, but due to sparse sampling and large distances between sites, such a result was very likely. It was suggested that R. tristis survived in an independent Pleistocene refugia both in the Alps and the Carpathians. Unfortunately, the area of the W Carpathians was not covered in this study and so the data from an important part of the species area of distribution was missing.
Western Carpathians represent the northernmost part of the Alpine-Carpathian mountain chain (Bielik, 1999). Due to high habitat diversity of this area, including numerous montane springs and streams, we hypothesize that the genetic patterns and level of the genetic diversity found by Bálint et al. (2011) at the large geographic scale could be also detected at a much smaller area of the W Carpathians. It could also be expected that this mountain range may include a contact zone between the Alpine (western) and Carpathian (eastern) lineages of the R. tristis. To test these assumptions, we generated and analysed a large dataset of R. tristis samples from the springs and streams in the W Carpathians and adjacent areas, which also allowed us to ll the gap in the phylogeographic puzzle of this species in the European mountains.

Results
The caddis y Rhyacophila tristis had a relatively wide distribution in the studied area of the Carpathians, both in streams and springs. The BIN algorithm and ASAP analysis resulted in identical partitioning of the 5' COI dataset into two MOTUs, whose distribution was largely disjunct (Fig. 2A).
The MOTU1, equal to BIN BOLD:AAD5574 (BIN_W), was recorded at 16 localities in the west, while the MOTU2, equal to BIN BOLD:ADL4166 (BIN_E), was found at 44 localities in the east. Both BINs were found only at two localities in the Fatra-Tatra area: the Biela Voda stream (T180) and Jazierce (V033). In the Fatra-Tatra area, where both BINs met, the highest levels of molecular divergence and a relatively clear boundary between the two BINS (MOTUs) were found ( Fig. 2A).
The comparison of the altitude range in the distribution of both BINs revealed signi cant differences. The BIN_W occurred at a signi cantly narrower altitude interval compared to BIN_E (p-value < 0.05; Fig. 2B).
On the other hand, no statistically signi cant preference for the habitat type (spring vs. stream) was detected for either BIN.
The haplotype network revealed the relationships between the two BINs identi ed in the study area as well as their connection with two additional R. tristis BINs from other European mountain ranges: the Balkan BIN_B (BOLD:ADL4367) occurring, e.g., in Rila Mts in Bulgaria and the BIN_A (BOLD:AAD5573) documented from the Swiss and Italian Alps (Fig. 2C). In total, 28 COI haplotypes of R. tristis were identi ed within 161 individuals from the W Carpathians. The BIN_W included 59 sequences grouped in 11 haplotypes, BIN_E 102 sequences in 17 haplotypes. There were found no signi cant differences in molecular genetic indices between the two BINs (Table 1).
Both BIN_W and BIN_E showed similar haplotype network patterns, i.e., star-like topology with a central, the most-frequent haplotype. RT01 was the most common haplotype in the BIN_W lineage, and was shared with German, Austrian, and Italian localities. The BIN_W included eight private haplotypes, six of them were found only in the Fatra-Tatra Area and two were private for the Western Beskids. Within the BIN_E, the total number of private haplotypes was 14, while four private haplotypes were identi ed in the Fatra-Tatra Area, the Slovak Ore Mts, and the Vihorlat Mts, one was found in the Poloniny Mts and the Apuseni Mts. The most common haplotype of the BIN_E (RT12 with 55 sequences) occurred at 32 sites.
Only a single haplotype (RT06) was shared among the Inner and Outer Eastern Carpathians (PM, VM, LB) and more distanced localities in Romania (Fig. 2C). Sequences from Romania were assigned to BIN_E, while sequences from Germany and Austria fell into the western BIN_W. The Bayesian phylogenetic reconstruction of R. tristis (5' COI) revealed that the divergence of the species likely started some 2.5 million years ago (Myr) and the western and eastern BINs diverged ca one million years later. Possible misidenti cation of sequences (marked as R. pubescens) downloaded from databases was also noted within the analysis (Fig. 2D). Both phylogenetic trees based on 3' or 5' COI sequences also showed that the BIN_E is more closely related to the Balkan BIN_B from the remote Rila Mts than to the BIN_W occurring in the same mountain system (Figs. 2D, S1).
The AMOVA of the whole dataset showed that most of the molecular variance could be attributed to the division of samples to geomorphological units or to river basins (Table S2). After partitioning the dataset into BIN_W and BIN_E, the highest proportion of the molecular variance was found within sampling localities (= subpopulations) in both BINs (Table 2). In the BIN_W, a relatively large portion of the total variance was explained by the differences among subpopulations within GU / RB. In the BIN_W, the variation at this level was negligible, however, in uence of the division into GU / RB was suggested ( Table  2).
The F ST values indicate a similar range of genetic differentiation within both BINs (Fig. 3). They showed that the most localities are relatively well connected in both BINs, but some level of isolation could not be refused among some subpopulations. The tests of isolation by distance revealed a positive correlation (Mantel test: r = 0.601; p-value = 0.001), supported also by the signi cant spatial autocorrelation (p-value = 0.000; Fig. S2). Such a positive and statistically signi cant correlation suggests a structuring effect of the geographical distance among studied localities.
The mismatch distribution analysis suggests a recent population expansion for both BIN_W and BIN_E, a scenario depicted by a typical unimodal distribution (Fig. 4A). In the BIN_W, this was also supported by the non-signi cant value of the Harpending's raggedness index (r). The non-signi cant values of the Harpending's raggedness index (r) was also found for BIN_E subpopulations of springs and streams separately (Fig. 5A). However, the mismatch distribution graph for BIN_E of streams has the character of a population expansion, as evidenced by neutrality tests that have statistically signi cant negative values (Table 3). In the case of BIN_E springs, all values of the neutrality tests were negative, which indicates a population expansion, but only the result of the Fu´s Fs test was signi cant (  Table 3).
The eBSP also indicates population expansion in both R. tristis BINs in the study area (Fig. 4B), which likely started before ca 50-80 ka in BIN_W and ca 25-30 ka in BIN_E. While the population of BIN_W grew gradually, the BIN_E population size increased rapidly in the last 20-25 ka (Fig. 4B). In the case of BIN_E, the expansion was found in populations of both springs and streams (Fig. 5B). The resulting eBSP plots show sudden population expansion and twice the present effective population size of BIN_E in streams in contrast to springs, where the trend of population expansion was gradual and not so prominent (Fig. 5B).

Discussion
Rhyacophila tristis is an upstream species (Céréghino, Cugny & Lavandier, 2002) requiring welloxygenated water with low organic matter levels (Robert & Curtean-Bănăduc, 2005). The results of the rst European phylogeographic study on this species by Bálint et al. (2011) showed strong genetic differences between its western and eastern populations, with lineages that probably survived independently in the Pleistocene refugia in the Alps and in the Carpathians. The presence of the western (BIN_W) and the eastern lineage (BIN_E) was also detected in a much smaller area of the W Carpathians. Moreover, based on the genetic landscape approach, we proved that W Carpathians form a contact zone between these BINs. Reconstruction of phylogeny suggests that the two BINs separated relatively long ago (~ 2 Myr, late Pleistocene), indicating that in fact they may represent two separate species. The analysis also showed that the eastern BIN_E is more closely related to the Balkan samples (BIN_B, Bulgarian Rila Mts) than to the BIN_W co-occuring in the same mountain system (W Carpathians). The climatic conditions in the mid and late Pleistocene impacted the geographical distribution of the species (Janeau & Aulagnier, 1997), which could also lead to intraspeci c differentiation or even speciation. When most of the mountains were covered by ice, some populations of species moved from mountains to adjacent lowland areas and some of them stayed in the isolated areas with suitable environmental conditions. The study of Asellus aquaticus (Isopoda) also con rms the survival of populations in the periglacial areas (Sworobowicz et al., 2020). However, to make a sound conclusion, additional studies including more molecular markers would be needed.
Concerning altitude, BIN_W occurred at a signi cantly narrower altitude interval compared to BIN_E, where the altitude of the localities ranged from approximately 200 to 1700 m a.s.l. Generally, elevation gradients shape the dispersal and population structure in mountain species due to isolation of high-elevation sites, local adaptation to extreme environments, or both (Hughes et al., 2008). Distribution of freshwater species can be also very often explained by altitude, but also ow type or different ecological demands (Illies & Botosaneanu, 1963). Possible relationship between the altitude range and the distribution scale of R. tristis BINs in our study area was indicated since the BIN_E was more widespread (44 localities) unlike the BIN_W, which was constrained to a smaller area (16 localities). However, it cannot be ruled out that BIN_W is extended further to the west, because two 3' COI sequences of Balint et al. (2009) from the Western Alps demonstrably belong to BIN_W (Fig. S1). If other samples analysed in the future con rm this assumption, it is likely that BIN_W could also have a wider altitudinal range.
Many freshwater species possess endemic genetic lineages or at least private haplotypes occuring in different European mountain systems, including the W Carpathians. For example, the caddis y Drusus discolor very likely persisted in the Tatra Mts in numerous refugia over multiple glacial cycles (Pauls, Lumbsch & Haase, 2006). Biodiversity richness of the W Carpathians is also supported by studies of the cold-adapted gammarids Gammarus balcanicus ( The AMOVA of the whole W Carpathian population of R. tristis revealed distinct variation among geomorphological units or river basins. This is however very likely due separate phylogeography, i.e., history and distribution of both detected BINs. The results of AMOVA on the level of individual BINs were more consistent with respect to GU or RB division, which suggests that considering BINs as separate evolutionary units is more natural. In both BINs, the highest portion of the variability was assigned to differences within subpopulations (localities), suggesting relatively recent colonization of the studied area. In the western BIN_W, however, reasonable part of the genetic variability was attributed to differences among subpopulations within geomorphological units / river basins, similarly as in Elmis aenea from springs of the same area (Bozáňová et al. 2020). On the other hand, AMOVA of BIN_E revealed no variability among localities within GU/RB, but slight differences between GU/RB were detected. This was similar to the more uvial ri e beetle Limnius perrisi from the same study or the black y Simulium degrangei (Jedlička et al., 2012) studied also in the W Carpathians. F ST values of both BINs agreed with the results of AMOVA. The differences between W Carpathian R. tristis lineages in the amount of genetic variation between localities within the GU/RB could indicate that the degree of local isolation is lower and / or the dispersal ability could be higher in BIN_E than in BIN_W. The eastern lineage also occurred in a signi cantly wider range of altitudes compared to western one. Altitude in uences a series of other factors, such as the proportion of suitable habitats or migration patterns, and indirectly it shapes species population genetics (Manel et al., 2003). Moreover, high-elevation sites are often reservoirs of speci c part of the species' genetic diversity (Hughes, 2007;Finn et al., 2011;Vuataz et al., 2016), but their vulnerability is high and the loss of isolated populations from higher elevations could result in a loss of important, locally adapted genotypes (Schiffers et al., 2013). Vice versa, adaptation potential and long-term survival of species or genetic lineages are dependent on su cient intraspeci c genetic diversity (Spielman, Brook & Frankham, 2004;Frankham, 2005), and from this perspective BIN_E could be at an advantage because it looks to be more tolerant to a wider scale of altitudes, and thus also environmental conditions. Differences between lineages (BINs) and, conversely, similar patterns of distribution of genetic variability with, e.g., ri e beetles from the same area but also with other insect species (Bunn & Hughes, 1997;Hughes et al., 1998) con rm the need to study population genetics in detail at the level of natural evolutionary units. It may allow us to better understand the evolution and phylogeography of individual taxa, or to look for general trends in the development of fauna of entire ecosystems.
Beside genetic variability distribution, the demography of the R. tristis BINs differs too in the W Carpathians. The BIN_W started to expand roughly around 40-50 thousand years ago (ka), while BIN_E around 25 ka during the Last Glacial Maximum (LGM). At that time, latitudinal temperature gradient existed across Europe, when winter soil temperatures were 10-20°C colder in Central and Northern Europe and around 2-4°C in Southern Europe than today (Barron & Pollard, 2002). Such conditions could suit the up-stream caddis y R. tristis, and so its glacial expansion could be expected (Bennet & Provan, 2008). The population of the W Carpathian BIN_E grew sharply compared to BIN_W, and its present effective population size is around ve times higher. We assume that this difference could be linked to the different distribution of the BINs in terms of altitude. BIN_E with a signi cantly wider altitude range probably found higher potential for dispersal in various mountain ranges within the W Carpathians. In addition, when comparing the eBSP of spring and stream habitats of BIN_E, a gradual, longer lasting increase in springs was observed. This could be explained by stable environmental conditions (and thus also longer time available for survival) in springs regardless of climatic changes (Minshall & Winger, 1968;Odum, 1971;Butler & Hobbs, 1982;Cushing & Wolf, 1984;Glazier & Gooch, 1987;Gooch & Glazier, 1991). Differences in eBSP and AMOVA could indicate that BIN_W survived in the study area for a longer period and that genetic drift occurred after the invasion of individual sites, causing genetic differences among populations within the GU/RB. In contrast, for BIN_E, the results suggest that it spread to the W Carpathians more abruptly and relatively rapidly after LGM. Subtle differences in the variability among GU/RB suggest that individual units/basins were colonised by a part of the originally homogeneous population, which brought speci c, although not signi cantly different compositions of haplotypes. However, the shorter time of occurrence in the studied area has not yet allowed BIN_E to create differences between local populations, as very likely happened within BIN_W.

Conclusion
Page 8/26 The analysis of the caddis y R. tristis of springs and streams in the Western Carpathians showed that these habitats form an important and unique part of the freshwater ecosystems producing or harbouring reasonable intraspeci c genetic diversity, and it can be assumed that these habitats play an equal role in other areas. The comparison with co-occurring species revealed very interesting results regarding the differences or similarities of their genetic diversity as well as history. This study did not have the ambition to address the taxonomy of the studied species, but it is clear that more attention should be paid to this, because, as genetic and geographical data have suggested, the possibility that the two identi ed lineages (BINs) represent separate cryptic species cannot be ruled out. The eastern lineage (BIN_E) most likely colonized the W Carpathians from Southeast Europe, probably at the end of the LGM, while BIN_W appears to originate in areas west of the W Carpathians but it survived and developed its population in the study area for a longer time. What is certain, however, is that the W Carpathians is the place where both lineages meet. What role do W Carpathian springs and streams play in survival and preservation of genetic diversity of R. tristis, whether the lineages hybridize in the contact zone or what are their speci c ecological requirements, will need to be further studied including additional samples from a wider area of distribution or other genetic markers.
Nevertheless, the data gained shifted our knowledge on the populations of the studied caddis y and suggested patterns that could be common to other freshwater species. Such information could be feasible for protection of aquatic ecosystems and preserving European freshwater biodiversity.

Study area
The study area includes the W Carpathians, one of the major geomorphological units of the Carpathian mountain range (Fig. 1A). It occupies approximately one-third of the total area of the Carpathians (Kondracki, 1989). Most of the W Carpathian ranges are of moderate altitude, ranging mostly from 500 to 1300 m a.s.l., only sparsely exceeding 1500 m a.s.l. (Kondracki, 1989). Geologically, the W Carpathians are diversi ed, especially in the inner part of the arc, with various kinds of bedrock (Mesozoic limestones and dolomites, Paleozoic granites and metamorphic rocks or Tertiary volcanic rocks). The outer W Carpathians are built almost solely from ysch (Grecula, 1997). During the Pleistocene cycles, the area of the W Carpathians remained mostly unglaciated. The continental ice sheet was only located in the northernmost foothills of the Polish part during its maximal extent in the Middle Pleistocene and the mountain glaciers covered valleys in the highest ranges (above 1700 m a.s.l.) (Lukniš, 1964). The mountain glaciers completely disappeared around 8.500 years ago (Lindner et al., 2003). The precipitation in the area is determined by differences in altitude and geomorphological relief. Hydrologically, the W Carpathian rivers have a rain-snow regime with oods in springs and summers.
The studied localities include 15 springs and 43 streams situated mainly in Slovakia, partially in Czechia and Poland, within the Inner and Outer Western Carpathian geomorphological units/subunits. These localities were supplemented by sampling sites in the Inner Eastern Carpathians (Vihorlat Mts -VM, Poloniny Mts -PM) and the Outer Eastern Carpathians (Low Beskids -LB) ( Fig. 1B; Table S1).

Sampling and morphological identi cation
The qualitative sampling of macrozoobenthos was performed in 2016-2017, within the framework of a broader hydrobiological research focused on karst springs and diversity of the W Carpathian streams. It was carried out by a multi-habitat kick-sampling technique (Frost, 1971), using a hydrobiological handnet with a mesh size of 0.5 mm. In the eld, specimens were xed in 96 % ethanol. In the laboratory, the invertebrates were sorted into higher taxonomic groups using stereomicroscope, pre xed with clean ethanol and stored in a freezer at − 25°C. Samples were supplemented by additional specimens of R. tristis from the collection of the Plant Science and Biodiversity Centre, Slovak Academy of Sciences (Bratislava). The individuals of R. tristis were morphologically identi ed using the determination keys Sedlák (1980)
These sequences were obtained from the BOLD (http://www.boldsystems.org). All sequences used in this study are included in the publicly available BOLD dataset DS-SKRHYTRI.
Two methods were used to test R. tristis samples for presence of the deeper lineages or speciesequivalents. The rst was BIN algorithm in BOLD (https://v4.boldsystems.org), where every uploaded sequence is compared to all records and assigned to an existing or a newly created Barcode Index Number (BIN) (Ratnasingham & Hebert, 2007). The BIN system clusters are unique and include molecular operational taxonomic units (MOTU) that closely correspond to the species. The second approach utilized the Assemble Species by Automatic Partitioning (ASAP), a novel and powerful method to build species partitions, based on pairwise genetic distances from single locus sequence alignments (i.e., barcode data sets) (Puillandre, Brouillet & Achaz, 2021).
The haplotype data les and the diversity indices were generated in DnaSP v5.10 (Librado & Rozas, 2009) based on 5' COI data. We calculated haplotype diversity (H), nucleotide diversity (π), number of polymorphic sites (S) and average number of nucleotide differences (K) per individual BINs of R. tristis species within the W Carpathians. Statistical comparison of the genetic indices between BINs was computed with the Wilcoxon signed rank test for paired data in R v4.0.2 (http://www.r-project.org). The same test was used to compare the altitude range between the individual BINs of the study area.
Haplotype networks were reconstructed using the median-joining method (MJN) in PopART v1.7 (Leigh & Bryant, 2015). The networks include sequences outside the W Carpathians to explain the phylogeographic relationships and haplotype distribution in the broader context.
The phylogeny was reconstructed based on 5' COI haplotypes using Bayesian approach in BEAST v2.5 (Bouckaert et al., 2019) and, separately, for the 3' COI marker also including published sequences Bálint et al., 2011). As an outgroup, the European congeners R. aquitanica and R. carpathica were included. The nucleotide substitution model was set through bModelTest (Bouckaerd & Drummond, 2017). The tree prior was set to Birth-Death and a strict molecular clock was used following the Path Sampling Selection. The strict molecular clock was calibrated with the standard mitochondrial rate for arthropod COI equal to 0.0115 substitutions/site/Myr (Brower, 1994). Two runs of Markov chain Monte Carlo (MCMC), each 20 million iterations long and sampled every 1,000 iterations, were performed for both fragments. Runs were examined using Tracer v1.7 (Rambaut et al., 2018), and all sampled parameters achieved a su cient sample size (ESS > 200). Tree les were combined using Log-Combiner v2.5.2 (Bouckaert et al., 2019), with the removal of the non-stationary 25 % burn-in phase. The maximum clade credibility chronogram was generated using TreeAnnotator v2.5.2 (Bouckaert et al., 2019).
The spatial diversity pattern in the studied area of Carpathians was illustrated with a genetic landscape approach generated with Alleles in Space (AIS; Miller, 2005). The genetic landscape visualizes the abrupt transitions between populations and groups of populations characterized by divergent haplotypes. First, using the AIS software, genetic distances between localities were calculated based on all 5' COI sequences and connected into a network based on the Delaunay Triangulation. The genetic distance values were set in the midpoints of each connection in the network. The raw genetic distances acquired were interpolated afterwards and the matrix of the 'elevation' values, with their respective latitude and longitude coordinates, was then imported into QGIS 3.18 (https://qgis.osgeo.org) software to produce a genetic divergence surface image using the inverse distance weighted algorithm. The resulting image was plotted onto a relief map to create a nal map in which the hypsometric tints (red-blue) re ect the genetic distance between population pairs.
To test if spatial distance is structuring the molecular diversity of 5' COI fragments, we run two types of isolation by distance tests: Mantel test (Mantel, 1967) and general spatial autocorrelation test using Alleles in Space (AIS; Miller, 2005). Both tests analyse correlation between spatial and molecular distance. To assess the signi cance, tests were run with 1,000 permutations.
The population structure of R. tristis in the W Carpathians based on 5' COI was characterized by the analysis of molecular variance (AMOVA) and xation indices (F ST ) using Arlequin v3.5 (Exco er & Lischer, 2010). The AMOVA was used to estimate whether the observed genetic diversity may be attributed to the geographical or river basin partitioning of populations in three levels: among geomorphological units (GU) / river basins (RB), among sampling sites within GU / RB and within sampling sites. We also performed AMOVA and F ST separately for individual BINs of the R. tristis, to nd out the differences between them. 138 sequences of R. tristis (36 localities) were included to calculate the F ST index, localities with 1 sequence were excluded. To test the signi cance of covariance components and xation indices, 1,000 permutations were performed.
Historical expansion patterns of the R. tristis in the studied area of Carpathians, based on 5' COI were examined in three different approaches. First, Tajima's D (Tajima, 1989), Fu's Fs (Fu, 1997) and Fu and Li's D (Fu & Li, 1993) neutrality tests with 10,000 permutations were calculated in DnaSP to test the selective neutrality and population stability. Secondly, nucleotide mismatch distribution was performed in Arlequin v3.5 (Exco er & Lischer, 2010) to detect the demographic and spatial dynamics of population expansion history in the W Carpathians. The t between the observed and expected distributions was evaluated by the test statistics of goodness-of-t, including the sum of squared deviation (SSD) and Harpending's raggedness index (r). Additionally, because part of the individuals was found in multiple springs and streams, the mentioned approaches were also performed separately for both habitats.    Table S1).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. SupplementaryMaterial.docx