DNA barcodes evidence the contact zone of eastern and western caddisfly lineages in the Western Carpathians

The region of the Western Carpathians is, among other aspects, very important for survival and diversity of European freshwater fauna due to the presence of a large number of (sub)mountain springs and streams. However, these ecologically and faunistically diversified habitats are still understudied in the context of genetic diversity and population structure of their inhabitants. This study focuses on genetic diversity and distribution patterns of the caddisfly Rhyacophila tristis, common and widespread representative of mountain freshwater fauna. Analysis of the COI mitochondrial marker revealed presence of the western and eastern lineages, with samples from both lineages being grouped in BOLD (Barcode of Life Data System) into separate BINs (Barcode Index Numbers). Our data indicates that eastern lineage (BIN_E) is more closely related to the Balkan populations than to co-occurring western lineage (BIN_W), and that the contact zone of the lineages passes through the W Carpathians. The study revealed phylogeographic and demographic differences between lineages, supporting hypothesis of their evolutionary independence and specific ecological preferences. The obtained genetic data of the R. tristis population from W Carpathians improved our knowledge about population genetics of this aquatic species and can contribute to understanding the state and evolution of biodiversity of freshwater ecosystems in Europe.


Results
The caddisfly R. tristis had a relatively wide distribution in the studied area, both in streams and springs, although not in high numbers. The BIN (Barcode Index Number) algorithm implemented in BOLD (Barcode of Life Data System) and ASAP (Assemble Species by Automatic Partitioning) analysis resulted in identical partitioning of the 5′ COI dataset into two genetic lineages, whose distribution was largely disjunct (Fig. 2A).
The western lineage, equal to BIN BOLD:AAD5574 (BIN_W), was recorded at 16 localities in the west, while the eastern lineage, equal to BIN BOLD:ADL4166 (BIN_E), was found at 44 localities in the east. Both lineages (BINs) were found in sympatry only at two localities in the Fatra-Tatra area: the Biela Voda stream (T180) and Jazierce spring (V033). In the Fatra-Tatra area, where both BINs met, the highest levels of molecular divergence and a relatively clear boundary between the two genetic lineages were found ( Fig. 2A). In the studied area, the BIN_W was recorded in a narrower elevation interval compared to BIN_E (Fig. 2B), data from the additional 11 western sites situated above 800 m a.s.l. (Supplementary Table S1) did not reveal presence of R. tristis in higher elevation.
The haplotype network revealed clear relationships between the two BINs identified 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 the Bulgarian Rila Mts and the BIN_A (BOLD:AAD5573) from the Swiss and Italian Alps (Fig. 3A). In total, 28 5′ COI haplotypes of R. tristis were identified within 161 individuals from the W Carpathians. The BIN_W included 59 sequences grouped in 11 haplotypes, BIN_E 102 sequences in 17 haplotypes. No significant differences were found 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 topologies with a central, the most-frequent haplotype. The most common haplotype in the BIN_W lineage (RT01) was shared with German, Austrian, and Italian localities. This BIN 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, four private haplotypes were identified 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 distant localities in Romania (Fig. 3A). 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 divergence of the species likely started some 2.5 million years ago (Myr) and the western and eastern BINs diverged ca. one and a half million years later. Possible misidentification of sequences (marked as R. pubescens) downloaded from public databases was also noted within the analysis (Fig. 3B). 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. 3B, S1).
The AMOVA of the whole dataset showed that most of the molecular variance could be attributed to the division of samples to physiographic 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 dots indicate outliers and the whiskers extend to the extreme values of the data, calculated as ± 1.5 × IQD from the median. Wilcoxon signed-rank test supported the significant differences in elevation range between two BINs (p-value < 0.05). www.nature.com/scientificreports/ the differences among subpopulations within physiographic units (PU) or river basins (RB). In the BIN_W, the variation at this level was negligible, however, influence of the division into PU/RB was suggested ( Table 2). The F ST values indicate a similar range of genetic differentiation within both BINs (Fig. S1). They showed that localities are relatively well connected in both BINs, but some level of isolation could not be refused among several subpopulations. The tests of isolation by distance revealed a positive correlation (Mantel test: r = 0.601; p-value = 0.001), supported also by the significant spatial autocorrelation (p-value = 0.000; Fig. S2). Such a positive and statistically significant correlation suggests a structuring effect of the geographical distance among studied localities. www.nature.com/scientificreports/ 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-significant value of the Harpending's raggedness index (r). The values of Tajima's D, Fu's Fs and Fu and Li's D neutrality tests for both BINs were negative and statistically significant (Table 3), which corresponds well with the Mismatch distribution test and supports population expansion hypothesis in both BINs (Fig. 4A). Overall, the results of the neutrality tests of the entire R. tristis population had negative values and the values of the Fu's Fs and Fu and Li's D tests were statistically significant ( Table 3).
The eBSP also indicates population growth in both R. tristis BINs in the study area ( Fig. 4B), which likely started between 40-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).

Discussion
Rhyacophila tristis is a common and widespread upstream species 31 requiring well-oxygenated water with low organic matter levels 28 . The results of the first European phylogeographic study on this species 29 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. When most of the mountains were covered by ice, populations of some 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 confirms such survival pattern in the periglacial areas 32 .
Our results show that the presence of the western (BIN_W) and the eastern (BIN_E) lineage was also detected in the W Carpathians themselves, and based on the genetic landscape approach, we proved that W Carpathians form a contact zone between these two lineages (BINs). Reconstruction of phylogeny suggests that the two BINs separated relatively long ago (~ 2.5 Myr, early Pleistocene), indicating (along with p-distance at the level of 2.08%) that in fact they could represent two separate species. However, to make a sound conclusion on their taxonomy, additional studies including more molecular markers as well as morphological feature analysis involving the study of adults would be needed. 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).
Many freshwater species contain endemic genetic lineages or at least private haplotypes occuring in different European mountain systems, including the W Carpathians. For example, the caddisfly Drusus discolor very likely www.nature.com/scientificreports/ persisted in the Tatra Mts in numerous refugia over multiple glacial cycles 23 . Biodiversity richness of the W Carpathians is also supported by studies of the cold-adapted gammarids Gammarus balcanicus 4,5 , G. fossarum 14 or G. jazdzewskii 33 . Although no endemic lineages of R. tristis were found in the W Carpathians, we have detected the presence of 21 private haplotypes within BIN_W and BIN_E, which have not yet been confirmed anywhere else. The AMOVA of the whole W Carpathian population of R. tristis revealed distinct variation among physiographic units and/or river basins. This is, however, very likely due to separate phylogeography, i.e., history and distribution of both detected BINs. On the contrary, the results of AMOVA on the level of individual BINs were more consistent with respect to PU 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, a reasonable part of the genetic variability was attributed to differences among subpopulations within PU/RB, similarly as in aquatic beetle Elmis aenea from springs of the same area 34 . On the other hand, AMOVA of BIN_E revealed no variability among localities within PU/RB, but slight differences between PU/RB were detected. This was similar to the more fluvial riffle beetle Limnius perrisi from the same study or the blackfly Simulium degrangei 35 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 on the level of genetic variation between localities within the PU/RB could indicate a lower degree of local isolation and/or higher dispersal ability of BIN_E than BIN_W.  www.nature.com/scientificreports/ The lineages (BINs) of R. tristis from the W Carpathians differed concerning elevation intervals in which they were found. Based on our data, we cannot unambiguously confirm whether this is a stable autecological characteristic of the lineages, because R. tristis has a significantly larger area of distribution than the studied area and to confirm or refute the possible difference, we would need more samples. On the other hand, in all sampled localities with an elevation above 800 m a.s.l. in the western part of the studied area, R. tristis was missing, which also suggests that the western lineage (BIN_W) could be characteristic for lower altitudes. Elevation influences a series of other factors, such as the proportion of suitable habitats or migration patterns, and indirectly it shapes species population genetics 1 . Moreover, high-elevation sites are often reservoirs of a specific part of the species' genetic diversity [36][37][38] , but their vulnerability is high and the loss of isolated populations from higher elevations could result in a loss of important, locally adapted genotypes 39 . Vice versa, adaptation potential and long-term survival of species or genetic lineages are dependent on sufficient intraspecific genetic diversity 40,41 , and from this perspective, i.e., more tolerant to different environmental conditions in different elevations, BIN_E could be at an advantage. On the contrary, if the BIN_W preference for lower elevation is confirmed, this lineage deserves more attention from the point of view of its protection due to higher anthropogenic pressure on lower-situated springs and streams.
By analysis of population genetics, we found the presence of differentiated genetic lineages (BINs) in R. tristis and, conversely, similar patterns of distribution of genetic variability as in beetles from the same area 34 , but also in other insect species 36,42 . These findings emphasize the need to study populations in detail at the level of natural evolutionary units. This allows us to better understand the evolution and phylogeography of individual taxa and to look for general trends in the development of the fauna of entire ecosystems.
Beside genetic variability distribution, also the demography of the R. tristis lineages detected in the W Carpathians differs. The BIN_W started to expand roughly around 40-80 thousand years ago (ka), while BIN_E around 25 ka during the Last Glacial Maximum (LGM). In 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 43 . Such conditions could suit the up-stream caddisfly R. tristis, and so its glacial expansion could be expected 44 . The population of the W Carpathian BIN_E grew sharply compared to BIN_W, and its present effective population size is around five times higher. We assume that this difference could be linked to the different distribution of the BINs. Differences in eBSP and AMOVA could indicate survival of the BIN_W in the study area for a longer period and occurrence of genetic drift after the colonisation of individual sites, causing genetic differences among populations within the PU/RB. In contrast, the results for BIN_E suggest that it spread to the W Carpathians more abruptly and relatively rapidly after LGM. Subtle differences in the variability among PU/RB suggest that individual units/basins were colonized by a part of the originally homogeneous population, which has brought with it specific, although not significantly different compositions of haplotypes. However, the relatively short time of occurrence in the studied area has not yet allowed BIN_E to fix more apparent differences between local populations, as very likely happened within BIN_W.

Conclusion
The detailed population study of the R. tristis caddisfly from the (sub)montane springs and streams in the W Carpathians showed that these habitats form an important and unique type of the freshwater ecosystems producing or harbouring substantial intraspecific genetic diversity. It can be assumed that these habitats play an equal role also 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 species of interest, but more attention should be paid to this issue in the future, because, as genetic and geographical data have suggested, the possibility that the two identified 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 and it survived and developed its population in the study area for a longer time. The Western Carpathians are the place where the two lineages meet. What role do W Carpathian springs and streams play in survival and preservation of R. tristis genetic diversity, whether the lineages hybridize in the contact zone or what are their specific ecological requirements, will need further study including samples from a wider area of distribution and additional genetic markers.

Materials and methods
Study area. The study area includes the W Carpathians, one of the major physiographic units of the Carpathian mountain range (Fig. 1A). It covers approximately one-third of the total area of the Carpathians 45 . Most of the W Carpathian ranges are of moderate elevation, ranging mostly from 500 to 1300 m a.s.l., only sparsely exceeding 1500 m a. s. l. 45 . Geologically, the W Carpathians are diversified, 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 flysch 46 . During the Pleistocene cycles, this area 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.) 47 . The mountain glaciers completely disappeared around 8.500 years ago 48 . The precipitation in the area is determined by differences in elevation and geomorphological relief. Hydrologically, the W Carpathian rivers have a rain-snow regime with floods 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 physiographic 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 S3). www.nature.com/scientificreports/ Sampling and morphological identification. 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 49 , using a hydrobiological hand-net with a mesh size of 0.5 mm. In the field, specimens were fixed in 96% ethanol. In the laboratory, the invertebrates were sorted into higher taxonomic groups using stereomicroscope, prefixed with clean ethanol and stored in a freezer at − 25 °C. Samples were supplemented by additional specimens of R. tristis from the ZooLab collection of the Plant Science and Biodiversity Centre, Slovak Academy of Sciences (Bratislava). The individuals of R. tristis were morphologically identified using the determination keys Sedlák from 1980 50  Data analysis. The obtained sequences were edited in SEQUENCHER v5.1 and aligned using the MUS-CLE algorithm 56 in MEGA v7 57 . In total, the 5′ COI dataset consisted of 194 sequences of R. tristis. It included 161 sequences from the W Carpathians (46 localities), the Poloniny Mts (3 loc.), the Vihorlat Mts (7 loc.) and Low Beskids (1 loc.) (Fig. 1). They were supplemented by 33 reference sequences of R. tristis from outside the W Carpathians (5-Austria, 6-Bulgaria, 5-Germany, 2-Italy, 11-Romania, 4-Switzerland), used for the reconstruction of haplotype networks and phylogenetic tree. These sequences were obtained from the BOLD (http:// www. bolds ystems. org).
Two methods were used to test R. tristis samples for presence of the deeper (genetic/evolutionary) lineages. The first was the BIN algorithm implemented in BOLD (http:// www. bolds ystems. org), where every uploaded sequence is compared to all records and assigned to an existing or a newly created Barcode Index Number (BIN) 58 . The BIN system clusters are unique and include sequences that presumably represent a single 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) 59 .
The haplotype data files and the diversity indices were generated in DnaSP v5.10 60 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-proje ct. org).
Haplotype networks were reconstructed using the median-joining method (MJN) in PopART v1.7 61 . 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 62 and, separately, for the 3′ COI marker also including published sequences 29,54 . As an outgroup, the European congeners R. aquitanica and R. carpathica were included. The nucleotide substitution model was set through bModelTest 63 . 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 64 . Two runs of Markov chain Monte Carlo (MCMC), each 20 million iterations long and sampled every 1000 iterations, were performed for both fragments. Runs were examined using Tracer v1.7 65 , and all sampled parameters achieved a sufficient sample size (ESS > 200). Tree files were combined using Log-Combiner v2.5.2 62 , with the removal of the non-stationary 25% burn-in phase. The maximum clade credibility chronogram was generated using TreeAnnotator v2.5.2 62 . The same methodology was used in case of the non-barcoding 3′ COI marker phylogeny reconstruction (Fig. 3).
The spatial diversity pattern in the studied area of Carpathians was illustrated with a genetic landscape approach generated with Alleles in Space (AIS) 66 . 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.16 (https:// qgis. org/) software to produce a genetic divergence surface image using the inverse distance weighted algorithm. The www.nature.com/scientificreports/ resulting image was plotted onto a relief map to create a final map in which the hypsometric tints (red-blue) reflect 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 67 and general spatial autocorrelation test using Alleles in Space (AIS) 66 . Both tests analyse correlation between spatial and molecular distance. To assess the significance, tests were run with 1000 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 fixation indices (F ST ) using Arlequin v3.5 68 . 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 physiographic units (PU)/river basins (RB), among sampling sites within PU/ RB and within sampling sites. We also performed AMOVA and F ST separately for individual BINs of the R. tristis, to find out the differences between them. 138 sequences of R. tristis (36 localities) were included to calculate the F ST index, localities with only 1 sequence were excluded. To test the significance of covariance components and fixation 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 69 , Fu's Fs 70 and Fu and Li's D 71 neutrality tests with 10,000 permutations were calculated in DnaSP v5.10 to test the selective neutrality and population stability. Secondly, nucleotide mismatch distribution was performed in Arlequin v3.5 68 to detect the demographic and spatial dynamics of population expansion history in the W Carpathians. The fit between the observed and expected distributions was evaluated by the test statistics of goodness-of-fit, including the sum of squared deviation (SSD) and Harpending's raggedness index (r).
Finally, the extended Bayesian skyline plot (eBSP), implemented in BEAST v2.6.2 software package 62 , was used to show the fluctuations of R. tristis demography in the W Carpathians over time. The model of substitution and molecular clock were set up identical as in the case of the phylogeny reconstruction. For comparison, two runs of Monte Carlo Markov Chains (MCMC) were performed, each 40 million iterations long and sampled every 10,000 iterations. The runs were examined in Tracer v1.7 and all the parameters reached the effective sampling size (ESS) above 200. The eBSP plots were produced using R v4.0.2 software (http:// www.r-proje ct. org). Both runs showed complementary results, thus only one is presented.