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 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 (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 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 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).
DNA extraction and PCR amplification
DNA was extracted from a total of 161 individuals, two legs of each individual, using the Chelex protocol (Casquet, Thebaud & Gillespie, 2012), followed by PCR amplification of ca. 650 bp-long barcoding fragment of the mitochondrial cytochrome c oxidase subunit I (5’ COI) using the primer pair LCO1490 - HCO2198 (Folmer et al., 1994). The PCR was performed in a total volume of 25 µl containing 5 µl of 5× DreamTaq™ Buffer, 1.5 µl of Mg+ 2 (25 mM), 0.5 µl of each primer (concentration 5 mM), 0.5 µl of dNTP Mix (20 mM), 0.125 µl (0.625 U) DreamTaq DNA Polymerase, 11.875 µl ultra-pure H2O and 5 µl of DNA template. The PCR cycling consisted of a 2-min initial denaturation at 94°C, followed by 40 cycles of 94°C (40 s) denaturation, 46°C (40 s) annealing and 72°C (1 min) extension and termination at 72°C (10 min) for a final extension. A 4 µl aliquot of the PCR products were visualized by GoldView (Solarbio) in electrophoresis on a 1 % agarose gel and GelLogic imaging equipment to check PCR product quality and length. The PCR products were purified with Exo-FastAP Thermo Scientific and were sent for sequencing to Macrogen Europe Inc., Amsterdam.
Additionally, the non-barcoding 3’ fragment of the COI was amplified using primers Jerry and Pat (Simon et al., 1994) from six selected individuals (see below for rationale) to create a phylogenetic tree together with 19 sequences (456 bp-long) of R. tristis and 15 sequences of congeneric species R. aquitanica McLachlan, 1879 (8) and R. carpathica Botoșaneanu, 1995 (7) from Bálint et al. (2009; 2011) (Fig. S1).
Data analysis
The obtained sequences were edited in SEQUENCHER v5.1 and aligned using the MUSCLE algorithm (Edgar, 2004) in MEGA v7 (Kumar, Stecher & Tamura, 2016). In total, the 5’ COI dataset consisted of 194 sequences of R. tristis. It included 161 sequences from the Western 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. 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 species-equivalents. The first 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 files 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., 2009; 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 sufficient sample size (ESS > 200). Tree files 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 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 (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 significance, 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 fixation indices (FST) using Arlequin v3.5 (Excoffier & 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 FST 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 FST index, localities with 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 (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 (Excoffier & Lischer, 2010) 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). Additionally, because part of the individuals was found in multiple springs and streams, the mentioned approaches were also performed separately for both habitats.
Finally, the extended Bayesian skyline plot (eBSP), implemented in BEAST v2.6.2 software package (Bouckaert et al., 2019), 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-project.org). Both runs showed complementary results, thus only one is presented.