Background {#Sec1}
==========

Biogeographic studies revealed wide ranges for many freshwater invertebrate species \[[@CR1]\]. This holds true in particular for species found in temperate and more northern latitudes, which had to recolonize habitats after glacial periods. Examples range from various insect species with terrestrial adult life stages to purely aquatic species, such as amphipod crustaceans and molluscs. Considering the typically wide ranges and high regional and local abundances, freshwater invertebrates are often regarded as frequent and long distance dispersers (e.g. \[[@CR2], [@CR3]\]). Over the past 15 years, however, the paradigm of wide ranges has been increasingly questioned, as molecular studies revealed the presence of morphologically cryptic species in many freshwater invertebrate taxa (e.g. \[[@CR4]--[@CR6]\]). These cryptic species often show rather small and allopatric ranges (e.g. \[[@CR7], [@CR8]\]), instead of the presumed broad distribution of the whole cryptic species complex.

One ecologically important freshwater taxon, which was formerly thought to be widely distributed in central and southeastern European streams, is the amphipod crustacean *Gammarus fossarum* KOCH, 1836. Several recent studies revealed an almost exponentially increasing number of overlooked species within the *G. fossarum* species complex with enhanced geographic sampling and improved sensitivity of molecular detection methods \[[@CR9]--[@CR11]\]. The highest species diversity, by far, within the *G. fossarum* species complex was found in the southeastern part of the range, where most of the newly discovered species were local endemics with narrow ranges \[[@CR11]\]. However, the four central and western European species, in particular clade 11, still show broad distributions \[[@CR10]\].

Generally, *G. fossarum* is mainly found in the upper reaches of streams and is sensitive to organic pollution \[[@CR12], [@CR13]\], high ammonium concentrations \[[@CR14]\], a lack of oxygen, and acidification \[[@CR15]\]. Owing to its high abundances and sensitivity to anthropogenic stressors, *G. fossarum* is often used in ecotoxicological studies (e.g. \[[@CR16]--[@CR18]\]). However, the precise cryptic species used in these experiments and whether a single or multiple species are used are rarely tested or reported. Validating species assignments prior to experiments is critically important, as studies explicitly investigating *G. fossarum* type A and B (here referred to as clade 12 and 11, after Weiss et al. \[[@CR10]\]) revealed ecological differences between the species \[[@CR19]--[@CR21]\]. In further studies comparing these two species, clade 11 was found to be more tolerant against tested stressors \[[@CR14], [@CR22]\], occurred in areas with higher human impact \[[@CR19]\] and was the better competitor in comparison with clade 12 \[[@CR13]\], but it also showed higher infection rates for various parasites \[[@CR23]\]. Additionally, in a direct comparison, populations of clade 11 were less differentiated across hundreds of kilometers than populations of clade 12 \[[@CR24]\], but still significant differentiation within clade 11 was found on a regional scale. These findings agree well with the moderate genetic differentiation found in a broad geographic area for members of clade 11 (e.g. \[[@CR10]\]). Even though these findings may indicate a relatively good dispersal ability for *G. fossarum* clade 11, it is difficult to predict actual dispersal rates, as they can be influenced by area- and site-specific environmental factors, like water chemistry, stream bed structure, land use and urbanization in the riverine environment, and fragmentation of streams by in-stream barriers, like dams or reservoirs (e.g. \[[@CR25]--[@CR27]\]). However, as understanding the patterns and mechanisms of dispersal and connectivity is crucial for predicting population resilience and long-term adaptability of a species \[[@CR28]\], it is important to determine the actual dispersal rates. An already regularly applied approach for this purpose is the use of genetic markers to estimate effective gene flow between populations, i.e. successful dispersal leading to genetic exchange between populations (e.g. \[[@CR3], [@CR29], [@CR30]\]).

In this study, we tested for factors driving the genetic structure of *G. fossarum* clade 11 in a human-impacted landscape at local and regional scales. To determine the population structure, we used two different genetic markers. For the main analyses, we used the barcoding fragment of the mitochondrial cytochrome c oxidase 1 (*CO1*) gene. We also examined nuclear microsatellite markers \[[@CR31], [@CR32]\] for a subset of populations to validate the *CO1* results. The study area was the Sauerland region, a low mountain range in North Rhine-Westphalia, Germany, which contains several small nature reserves, but is also used for agriculture, industry, forestry, and tourism. The hydrological structure of streams in the Sauerland region is strongly influenced by anthropogenic factors, such as in-stream barriers occurring approximately every 1,000 m \[[@CR33]\]. Therefore, the region is characterized by high site heterogeneity in terms of ecological parameters as well as habitat fragmentation, making it an interesting area to study the impact of anthropogenic factors on the realized dispersal of aquatic invertebrates. To account for these factors, we characterized sampling sites based on several ecological parameters and combined dense small-scale sampling with broader regional sampling within a range of 85 km. Specifically, we tested the following hypotheses:Populations of *Gammarus fossarum* clade 11 are genetically differentiated at the regional scale when considering the whole study area. In contrast differentiation is low at the local scale of few kilometers with evidence of gene flow between populations. These expectations are based on previous genetic analyses of *G. fossarum* species (e.g. \[[@CR24], [@CR34]--[@CR36]\]).Genetic variation is significantly partitioned after drainage basins at any spatial scale within the riverine network according to the Stream Hierarchy Model \[[@CR37], [@CR38]\], because *G. fossarum* is confined to aquatic habitats throughout its life cycle. Further, populations show patterns of isolation by distance (IBD) within catchments, especially when considering the waterway distance as this pattern was also found in other population genetic studies of *G. fossarum* clade 11 (e.g. \[[@CR24]\]).

After testing the two hypotheses, we discuss the findings in order to identify possible drivers of population structure.

Methods {#Sec2}
=======

Study site and sampling {#Sec3}
-----------------------

Specimens of *G. fossarum* were collected via kick-sampling at 54 sampling sites in 45 different streams in the Sauerland region (Germany). Animals were preserved in 96 % ethanol. The sampling sites were located in three major catchments (Ruhr: 48; Eder: 5; Lippe: 1). For the analyses, the Ruhr catchment was further divided into four sub-catchments, Lenne (12), Möhne (5), Volme (2), and Ruhr (29), where Ruhr means that the streams flow directly into the Ruhr and not via another main river, such as the Lenne (see Table [1](#Tab1){ref-type="table"} for details). To characterize the sampling sites, various parameters were measured (see Additional file [1](#MOESM1){ref-type="media"}). Directly determinable parameters were coordinates, sub-catchment association, altitude, distance to the spring, and if the sampling site was positioned in a small nature reserve. Land use types (estimated as percentages) in a buffer zone of 10 m wide and 1 km long upstream were determined in QGIS v. 2.8.2 \[[@CR39]\] using the ATKIS land cover vector data \[[@CR40]\]. The categories were conifer, broadleaf, and mixed forests, farmland, grassland, water bodies, and urban areas. The hydromorphological variables channel pattern (Main Parameter 1, MP1), longitudinal profile (MP2), and structure of water beds (MP3) were derived from the national hydromorphological survey (\[[@CR41]\], described in \[[@CR42]\]). Further, the ecological status according to the Water Framework directive \[[@CR43]\] was determined for the sampled streams. Data for both hydromorphological variables and ecological status were provided by the federal state authority LANUV (Landesamt für Natur, Umwelt und Verbraucherschutz © Land NRW, Recklinghausen, [http://www.lanuv.nrw.de](http://www.lanuv.nrw.de/)). Chemical measurements were obtained from the ELWAS-web portal. Most of the official measuring points were not directly located at the sampling sites in this study, but 41 of these measuring points were located within a 5 km reach. The frequency with which the chemical values were measured differed among sites, varying from one to four years and one to seven times a year, but not all parameters were measured on all dates. To correct for these differences and variation within sampling sites, the average value for each chemical was calculated. The following chemical parameters were used for the analyses: calcium, iron, oxygen, chloride, ammonium, total organic carbon (TOC), total nitrogen, and pH (for details see Additional file [1](#MOESM1){ref-type="media"}).Table 1Sampling sites and number of studied *G. fossarum* specimens (n) in the Sauerland areasitestream namenlongitude (GK3)latitude (GK3)sub-catchmentgenetic markersampling yearHSK10Hannebecke123458921.35685865.4RuhrCO12014HSK6bNierbach263455373.85686737.7RuhrCO12014KL14Bieberbach603426167.25696528.9RuhrCO12014KL3Waldbach713432816.05681103.7RuhrCO12014KL9Röhr 3593434350.55682378.1RuhrCO12014LELeiße533452683.75673900.3RuhrCO1 & msat2011 & 2013MEMedebach183466979.75687670.3RuhrCO12012NGRenau463462376.45675757.3RuhrCO1 & msat2011 & 2013NLNamenlose123464929.25675391.8RuhrCO12011PL1Palme 1253457449.95680115.8RuhrCO1 & msat2011PL2Palme 2823458169.15676718.3RuhrCO12011 & 2013QB12Elpe363461411.55681943.7RuhrCO12013QB17Ilpe 1363445685.85678145.5RuhrCO12013QB22Kleine Henne583453172.35688126.2RuhrCO12014QB23Ilpe 2363447494.95677477.6RuhrCO12014QB24Hengsbecker Bach483442318.15677813.1RuhrCO12014QB27Schürenbach583446183.35688992.2RuhrCO12014QB29Krähe463425599.15681574.1RuhrCO12014RO1Röhr 1203434317.55679605.7RuhrCO12012RO2Röhr 2253431365.35686986.9RuhrCO12012RU3Ruhr 3523466754.45681626.1RuhrCO1 & msat2011 & 2013RU4Ruhr 473467531.05676847.1RuhrCO12011SBSchlebornbach583460025.75694421.7RuhrCO1 & msat2011 & 2013VA1Valme 1203459155.35678350.2RuhrCO12011VR11Refflingser Bach473406403.25698078.2RuhrCO12013VR16Elsebach363404615.05697289.2RuhrCO12013VR17Palme 3443457758.05678654.0RuhrCO12014VR5Valme 2483458707.35683774.6RuhrCO12013VR7Kelbke483445543.35686405.8RuhrCO12013E01Mühlenbach143410770.95681556.6Lenne (Ruhr)CO12011E02Gleierbach543453326.75672099.4Lenne (Ruhr)CO1 & msat2011 & 2013E04Husberger Bach293409316.45681432.9Lenne (Ruhr)CO12011E06Fretterbach153434941.85674698.5Lenne (Ruhr)CO1 & msat2011E11Elspe223438640.35672688.0Lenne (Ruhr)CO12011GBGrüner Bach293408046.75691018.5Lenne (Ruhr)CO1 & msat2012KL13Schwarze Ahe683410745.15675226.3Lenne (Ruhr)CO12014KL15Krummenau593409854.35660639.6Lenne (Ruhr)CO12014KL2Worbscheider Bach573417114.15663534.7Lenne (Ruhr)CO12014NBNimmer Bach373400329.55688143.1Lenne (Ruhr)CO12011 & 2012SOSorpe 1143458195.65673861.6Lenne (Ruhr)CO12011VR2Sorpe 2563460431.05673996.4Lenne (Ruhr)CO12013AAAa173467971.25696423.1Möhne (Ruhr)CO12011GSGroße Schmalenau303441171.55701881.2Möhne (Ruhr)CO12012KL6N.N.603465023.55697739.9Möhne (Ruhr)CO12014LOLörmecke303458555.95702556.6Möhne (Ruhr)CO12011QB10Hirschberger Bach273453285.95701355.8Möhne (Ruhr)CO12013VR12Ennepe483395008.05671675.1Volme (Ruhr)CO12013VR23Epscheider Bach563393785.45683027.4Volme (Ruhr)CO12014BBBremke-Bach293466644.65668828.9EderCO12011HB1Hallebach 1163478675.05677590.1EderCO12011HB2Hallebach 2203474568.95675458.9EderCO12011NDNeerdar243478218.45682425.6EderCO12011NHNuhne583467984.55671089.2EderCO1 & msat2011 & 2013ALAlme303473839.95701645.1LippeCO1 & msat2011

*G. fossarum* sampling was conducted in 2011, 2012, 2013, and 2014 (Table [1](#Tab1){ref-type="table"}). Samples from six sites in 2011 were provided by Maria Gies and Martin Sondermann (University Duisburg-Essen). To examine changes in the genetic structure of populations over time, seven 2011 sampling sites were resampled in 2013 and one 2011 site was sampled again in 2012 (Table [1](#Tab1){ref-type="table"}). Additional samples from 2013 and 2014 were originally obtained for another study in which a stream was sampled at four sites every 200 m with an in-stream barrier between separating sites two and three. At seven sampling sites in 2014, a similar sampling scheme was used for sewage plants, and, at two sampling sites, similar schemes were used for mining sites. Since no differences in haplotype frequencies were found between those sites (*F*~ST~ values were not significant, data not shown), the samples were merged into a single sampling site for this study.

Sequencing and genotyping {#Sec4}
-------------------------

DNA was extracted from the pereopods of 2,086 specimens using a modified salt-extraction protocol \[[@CR44]\] (see Additional file [2](#MOESM2){ref-type="media"} for details of used protocol). For all specimens, a fragment of the mitochondrial barcoding gene *CO1* was amplified with the standard primers HCO2198 and LCO1490 \[[@CR45]\] using the following PCR protocol: 1× PCR buffer, 0.2 mM dNTPs, 0.5 μM each primer, 0.5 μl of DNA template, 0.025 U/μl HotMaster Taq (5 PRIME GmbH, Hamburg, Germany), filled to 15 μl with sterile H~2~O. The PCR settings for *CO1* amplification were as follows: initial denaturation at 94 °C for 2 min; 34 cycles of denaturation at 94 °C for 20 s, annealing at 46 °C for 30 s, and extension at 65 °C for 60 s; final extension at 65 °C for 5 min. The PCR products (9 μl) were purified using 0.5 μl of ExoI (20 U/μl) and 1 μl of FastAP (1 U/μl, both enzymes; Thermo Fisher Scientific, Schwerte, Germany). The reaction was incubated for 25 min at 37 °C followed by an inactivation step at 85 °C for 15 min. Purified PCR products were bidirectionally sequenced by GATC-Biotech (Konstanz, Germany).

For 420 specimens from ten sampling sites (Table [1](#Tab1){ref-type="table"}), eight microsatellite loci were additionally amplified: Gam 2, Gam 14 \[[@CR31]\], Gamfos 10, Gamfos 13, Gamfos 18, Gamfos 22, Gamfos 27, and Gamfos 28 \[[@CR32]\]. For each primer pair, the originally published sequence was adapted to our analysis system by adding a universal M13 tail (5′-CAC GAC GTT CTA AAA-3′) to the 5′ ends of the forward primers for primers developed by Danancher et al. \[[@CR31]\] and to the reverse primers for those developed by Westram et al. \[[@CR32]\]. The amplification of Gam 2, Gam 14, Gamfos 13, Gamfos 27, and Gamfos 28 was performed using the following protocol: 1× PCR buffer, 0.2 mM dNTPs, 0.2 μM sequence-specific untailed primer, 0.05 μM sequence-specific tailed primer, 0.2 μM fluorescently labelled universal M13 primer, 5 % dimethyl sulfoxide (DMSO), 0.5--1 μl of DNA template, 0.025 U/μl HotMaster Taq (5 PRIME GmbH, Hamburg, Germany), filled to 15 μl with sterile H~2~O. For Gamfos 10, 18, and 22, betaine was used instead of DMSO. PCR settings for the amplification of Gam2, Gam14, Gamfos 13, Gamfos 18, Gamfos 22, and Gamfos 28 were as follows: initial denaturation at 94 °C for 2 min; 36 cycles of denaturation at 94 °C for 20 s, annealing at 60 °C for 30 s, and extension at 65 °C for 30 s; final extension at 65 °C for 45 min. For Gamfos 10 and Gamfos 27, an annealing temperature of 54 °C was used.

Allele sizes were determined using polyacrylamide gel electrophoresis on a LI-COR® 4300 DNA Analyzer with the software Saga^GT^ (LI-COR Biosciences, Bad Homburg, Germany). For Gam 2, the alleles could not be scored because too many stutter bands occurred over a range of 90 base pairs.

Sequence data analysis {#Sec5}
----------------------

Sequences were assembled and edited in Geneious 8.1.2 ([http://www.geneious.com](http://www.geneious.com/), \[[@CR46]\]), and a multiple sequence alignment was computed with MAFFT \[[@CR47]\] as implemented in Geneious (automatic algorithm selection, default settings). Genetic diversity was calculated as haplotype diversity using Arlequin 3.5 \[[@CR48]\]. A minimum spanning network was calculated using Arlequin and visualized using HapStar 0.7 \[[@CR49]\]. To show the position of the haplotypes in a broader phylogenetic context nine sequences from five different clades (3, 10, 11, 12 and 13), which have been used in the study of Weiss et al. \[[@CR10]\] (GenBank Accession numbers: JN900490, KF521805, KF521817, KF521822, KF521828, KF521829, KF521832, KF521833, KF521835), were added to a dataset of the main haplotypes of this study and a neighbor-joining Tree was calculated with MEGA6 \[[@CR50]\], with evolutionary distances computed using the Kimura 2-parameter method. To test for population differentiation, pairwise *F*~ST~ values between populations from different sampling sites were calculated using Arlequin. Negative values were set to zero. The significance levels for this and all subsequent *F*~ST~ calculations were adjusted to account for multiple testing using the false discovery rate control (FDR, \[[@CR51]\]). Additionally, *F*~ST~ values were calculated using data from multiple years for populations originating from the same stream to test if haplotype frequencies changed significantly over time.

Microsatellite data analysis {#Sec6}
----------------------------

Using the program MICRO-CHECKER 2.2.3 \[[@CR52]\], our data set was checked for the occurrence of allelic drop out, stutter bands, and null alleles and null allele frequencies were calculated using the Dempster method \[[@CR53]\]. All loci in all populations were tested for deviations from Hardy--Weinberg equilibrium (HWE) and from linkage equilibrium using Arlequin. To estimate genetic diversity, allelic richness was calculated with the program HP-Rare 1.0 \[[@CR54]\] using rarefaction to correct for differences in sample sizes.

Genetic distances between populations were calculated as pairwise *F*~ST~ values in Arlequin. To correct for null alleles, FreeNA \[[@CR55]\] was used, in which the ENA correction was implemented. To test if the allele frequencies in the populations changed over time, *F*~ST~ values were also calculated between populations originating from the same stream, but sampled in different years.

Determination of drivers of genetic differentiation (*CO1* data) {#Sec7}
----------------------------------------------------------------

Different landscape genetic approaches were used to identify drivers of genetic differentiation in the study area. First, to test if the distribution of genetic variance can be explained by the partitioning of populations into the six sub-catchments, AMOVA (analysis of molecular variance, \[[@CR56]\]) was implemented in Arlequin. To test if genetic distances were correlated with geographic distances (isolation by distance, IBD), Mantel tests were performed in Arlequin. Pairwise *F*~ST~ values were used as a measure of genetic distance, and either straight-line or waterway distances were used for geographic distance. Straight-line distances were calculated for all population pairs and waterway distances for populations within the Ruhr catchment (*n* = 48) using QGIS v. 2.8.2 \[[@CR39]\] and the map containing the streams provided by the federal state authority LANUV (Gewässerstationierungskarte des Landes NRW © LANUV NRW (2013)). To test if in-stream barriers influenced the connectivity between populations, another Mantel test was conducted with 21 populations of the Ruhr sub-catchment using pairwise *F*~ST~ values and the number of in-stream barriers between population pairs per km waterway distance, first for barriers \>0.5 m and then for barriers \>1 m. The number of in-stream barriers was calculated from the QUIS database \[[@CR57]\]. Barriers were mostly weirs and barrages, often classified as only partly or not at all passable for the aquatic fauna in the ELWAS-web portal. To investigate if spatial factors other than geographic distance or catchment assignment shape genetic structure and to infer the number of populations, spatial Bayesian clustering models were implemented in the R package GENELAND, v4.0.5 (\[[@CR58], [@CR59]\], R version 3.2.2, \[[@CR60]\]). The number of populations K was allowed to vary between 1 and 20. For six independent runs, 1,000,000 MCMC iterations were calculated, sampling every 100 steps. The maximum number of nuclei in the Poisson-Voronoi tessellation was set to 6,300. For post-processing, a burn-in of 1,500 iterations was used and the pixels along the *X*-axis were set to 300 and along the *Y*-axis to 150 according to the ratio of the sampling area. A second AMOVA implemented in Arlequin using the groups detected in GENELAND was used to analyze if this clustering better reflects the population structure than the sub-catchments.

To test which factors determine the association of the individual populations to the GENELAND groups, a discriminant analysis was performed in IBM SPSS v. 23. For this analysis, the environmental variables obtained to characterize sampling sites were independent variables and the main GENELAND groups (consisting of \>3 sites) were dependent variables. Additional to the environmental variables, also the geographic position of the sampling sites, represented by their coordinates, was used as a variable to include a proxy for the possible influence of colonization history on genetic structure. To standardize the predictor variables, z-scores were calculated in SPSS prior to the discriminant analysis. Further, multicollinearity between the 25 environmental variables was tested by calculating the variance inflation factors (VIFs) for each variable using the R package "usdm" \[[@CR61]\]. Collinearity between variables was assumed for VIFs \>10. After calculating the VIFs, the variables, which were highly collinear, were stepwise excluded from the analysis until all values were \> 10. Therefore, the following discriminant analysis was conducted with a reduced set of variables. In the discriminant analysis, the stepwise Wilks lambda method was used to select predictors.

Results {#Sec8}
=======

We generated *CO1* sequences for 2,086 *G. fossarum* clade 11 specimens from 54 sampling sites. In the 557 bp alignment, we detected 52 variable sites, of which 12 were non-synonymous substitutions. Specimens clustered into 40 distinct haplotypes (H1--H40), of which ten had a frequency of over 1 % (H1, H7, H12, H19, H21, H28--31, and H33) and were found in 97.2 % of all specimens (Additional file [3](#MOESM3){ref-type="media"}). We observed the most common haplotype, H1, in 38 populations and 48.8 % of sequences. Other common haplotypes were H19 (19 populations; 16.7 % of sequences), H21 (8; 6.8 %), and H33 (6; 11.6 %), and we found all four haplotypes in at least three sub-catchments. Each of these common haplotypes had further derived ones, which were mostly differentiated by a single mutation, visible in the haplotype network (Fig. [1](#Fig1){ref-type="fig"}). H33 and the surrounding haplotypes H34 to H40 were differentiated by at least 14 substitutions (distance between H33 and H31) from all other haplotypes. In the neighbor-joining tree (Additional file [4](#MOESM4){ref-type="media"}), all haplotypes clustered together with sequences belonging to Clade 11. As already visible in the network, this clade was divided into two sub-clusters (one containing H33 and the second all other main haplotypes), but the distance between the two sub-clusters was shallow as compared to distances between the different clades.Fig. 1Minimum spanning network created from *CO1* sequences. Circles represent different haplotypes and their dimensions are scaled based on the number of sequences, which are given in Table [1](#Tab1){ref-type="table"}. Vertical lines represent missing or unsampled haplotypes. Red edges of circles indicate that these haplotypes were found at different sampling sites, while black edges indicate private haplotypes. Haplotypes are colored similar to Fig. [2a](#Fig2){ref-type="fig"}

We only detected eight haplotypes that were shared between at least two populations (H1, H7, H19--22, H31, and H33) and we observed private haplotypes in 20 populations. At two sampling sites, i.e., KL15 (*n* = 59) and KL2 (*n* = 57), populations consisted exclusively of one private haplotype each. Overall, we detected between one and six haplotypes per population, with an average of two haplotypes, and a haplotype diversity of between 0.00 and 0.68, with an average of 0.16 (Additional file [3](#MOESM3){ref-type="media"}). The nucleotide diversity ranged from 0.000 to 0.136 with an average of 0.048 (Additional file [3](#MOESM3){ref-type="media"}). In seven of the eight populations for which we compared haplotype frequencies among years, we did not observe significant differences. The only exception was population NH, where, in the second sampling year, only the main haplotype of the population (H21) was rediscovered together with a new haplotype (H20), resulting in a small, but significant, *F*~ST~ value.

Of the seven microsatellites analyzed to complement the *CO1* data set, only four were polymorphic in the studied populations (Gamfos 10, 13, 18, and 28). Three of the ten populations were also monomorphic for the same allele at Gamfos 13. We observed between 4 and 23 alleles in all populations at the different loci. We did not detect significant linkage between loci considering all populations. We found evidence for null alleles for three of the loci in different populations. The null allele frequencies, number of alleles, and observed and expected heterozygosity (H~O~ and H~E~) are given in Additional file [5](#MOESM5){ref-type="media"}. Null allele frequencies ranged from 0.00 to 0.29. We observed deviations from HWE in all populations for at least one locus (see Additional file [5](#MOESM5){ref-type="media"}). As a measure of genetic diversity, we estimated allelic richness using rarefaction for a minimum sample size of 15 diploid individuals. We observed an average allelic richness over all loci of between 2.5 and 5.3 and we detected private alleles in six of the populations (Additional file [5](#MOESM5){ref-type="media"}). In two of the populations sampled twice in different years (NH and E02), the allele frequencies changed over time, resulting in small, but significant, *F*~ST~ values.

Regional and local differentiation {#Sec9}
----------------------------------

The *CO1* haplotype composition at the sampling sites differed strongly in many cases on a regional as well as on a local scale (Fig. [2a](#Fig2){ref-type="fig"}), resulting in an overall high differentiation (mean *F*~ST~ = 0.61), and 81 % of all pairwise *F*~ST~ values were significant (1,152 out of 1,431 comparisons; Fig. [2b](#Fig2){ref-type="fig"} and Additional file [6](#MOESM6){ref-type="media"}). Considering only the subset of populations in which microsatellites were analyzed, 82 % of the *F*~ST~ values indicated significant population differentiation. All of the populations in this subset were significantly differentiated from each other when analyzing microsatellite loci (Fig. [2b](#Fig2){ref-type="fig"}, Additional file [7](#MOESM7){ref-type="media"}), with *F*~ST~ values ranging between 0.07 and 0.57. When estimating *F*~ST~ values with FREENA using the ENA method to correct for null alleles, they were slightly lower in most cases, ranging from 0.06 to 0.55 (Additional file [7](#MOESM7){ref-type="media"}). Altogether, we detected strong signatures of local and regional isolation of *G. fossarum* populations using both data sets*.*Fig. 2**a** *CO1* haplotype map showing the haplotype composition for *G. fossarum* at different sampling sites in the Sauerland region. The sizes of haplotype pie charts are scaled according to the numbers of sequences per site, which are given in Table [1](#Tab1){ref-type="table"} together with the sub-catchment association of sampling sites. Red stars indicate water reservoirs. Red highlighted sampling sites indicate that microsatellites were analyzed at these sites. All private low-frequency haplotypes are colored in black, or in gray if more than one private haplotype was found at the respective site. Colored contour lines illustrate GENELAND groups, named A--H. **b** Bar chart showing the frequency of significant and non-significant *F* ~ST~ values for microsatellites (msat) and *CO1* sequences

Drivers of local isolation {#Sec10}
--------------------------

The AMOVA using six sub-catchments as groups revealed significant population differentiation between sub-catchments, with 19.6 % of the variation partitioned between groups. However, 59.8 % of the variation was between populations within sub-catchments (Table [2](#Tab2){ref-type="table"}). Mantel tests revealed small but significant correlations between genetic distance (pairwise *F*~ST~) and geographic distance. Here, the fit was better using straight-line distances (*R*^2^ = 0.17, *p* = 0.000) than waterway distances (*R*^2^ = 0.03, *p* = 0.010, see Fig. [3](#Fig3){ref-type="fig"}). We did not detect correlations between genetic distance and the number of barriers per km between a subset of 21 sampling sites, for two barrier size thresholds (barriers \>0.5 m: *R*^2^ = 0.0003, *p* = 0.433; barriers \>1 m: *R*^2^ = 0.0007, *p* = 0.408).Table 2Results of AMOVA analysis according to 1) sub-catchments and 2) GENELAND groupsBetween sub-catchmentsBetween GENELAND groupsSource of variationd.f.% variationFixation indexd.f.% variationFixation indexbetween groups*F* ~*CT*~519.6**0.20**672.1**0.72**between populations within groups*F* ~*SC*~4859.8**0.74**4710.0**0.36**within populations*F* ~*ST*~203220.6**0.79**203217.9**0.82**d.f. = degrees of freedom; bold values for Fixation index indicate significant population differentiationFig. 3**a** Correlation between pairwise genetic and waterway distances for sampling sites of the whole Ruhr catchment. **b** Correlation between pairwise genetic and straight-line distances for all sampling sites

Other major barriers located between some of the sampling sites were water reservoirs (see Fig. [2a](#Fig2){ref-type="fig"}). One such water reservoir, the Bigge reservoir, separated populations KL2 and KL15, and a second smaller one (Ahauser reservoir) separated both sites from all other sampling sites; both populations consisted only of one private haplotype each (H28 and H29, resp., see Fig. [2a](#Fig2){ref-type="fig"}). Another reservoir (Sorpe reservoir) separated QB29 from the other sampling sites. Between QB29 and the next sampled population (RO2), two haplotypes were shared, but the haplotype composition was significantly different. The population GS of the Möhne sub-catchment was separated by the Möhne reservoir from the other populations of this catchment. This reservoir also separated the whole sub-catchment from the Ruhr catchment. In both the GS population and some populations from the Ruhr and Möhne catchment (e.g. KL9, VR7, and LO), we only detected haplotype H1, and we did not detect a barrier effect of the Möhne reservoir based on the *CO1* sequence analyses.

To infer the population structure more directly, without manually assigned assumptions regarding affiliations to catchment areas, we performed a clustering analysis with GENELAND. In this analysis, we observed seven distinct groups, termed A--G (letters in circles in Fig. [2a](#Fig2){ref-type="fig"}). The biggest group was group A, with 24 populations containing mostly or exclusively haplotype H1. The second biggest group, dominated by H19, was group B with 13 populations. In most of these populations, we also detected H1, although at minor frequencies. This group was geographically split into two subgroups, with two populations located in the western area, and most populations in the east, with no visible connection between the two subgroups. Groups C and D were the easternmost and westernmost groups, both consisting of six populations, dominated by H21 and H33, respectively. While populations of group C shared haplotypes with groups A and B, all haplotypes of group D (H33 to H40) were exclusively found in that group. All but one population of group D contained additional private haplotypes, aside from the main haplotype. Group E consisted of three populations that shared H1 at a minor frequency, but had different dominant haplotypes, which were all derived from H1 (Fig. [1](#Fig1){ref-type="fig"}). Groups A to E all contained populations from at least two sub-catchments. The last two groups F and G each consisted only of one population (KL2 and KL15) and contained exclusively one private haplotype (H28 and H29, respectively). Both haplotypes were three mutations apart from H19 and six from each other. Using the GENELAND groups in an AMOVA, we found that this clustering reflects the population structure better than the sub-catchment structuring, as 71.1 % (instead of 19.6 %) of the variation was between groups and only 10 % (instead of 59.8 %) was among populations within groups (Table [2](#Tab2){ref-type="table"}).

To test if other environmental factors determine the observed population structure, we gathered data for 25 environmental variables (Additional file [1](#MOESM1){ref-type="media"}). As groups F and G only consisted of single sites and populations of group E were highly differentiated from each other, we only conducted the following analyses for the four main groups, A--D. Testing for multicollinearity revealed that some of the variables were highly correlated. We therefore excluded the following variables for subsequent analyses: altitude, MP3, grassland, total nitrogen, and pH, resulting in 20 variables for the discriminant analysis. The stepwise discriminant analyses revealed that the most useful variable for predicting assignments to the GENELAND groups was longitude; all other variables were not included in the model. With the resulting model, 65 % of the populations clustered correctly to the four groups, with differences in prediction performance between groups (see Additional file [8](#MOESM8){ref-type="media"}). While all populations of groups C and D clustered correctly, this was the case for only 61.5 % of group B and 50.0 % of group A populations, as these two groups had a broader geographic range with overlapping longitudinal values.

Discussion {#Sec11}
==========

We analyzed a large number of specimens and populations of the freshwater crustacean *G. fossarum* clade 11 in the anthropogenically heavily impacted Sauerland region to identify factors that influence population structure and limit dispersal in this species.

Our first hypothesis was that populations are genetically differentiated when considering the regional scale of the whole study area, whereas the differentiation is low at the scale of few kilometers. In agreement with the first expectation of this hypothesis we detected strong regional differentiation in the Sauerland area, especially between the eastern and westernmost populations. The westernmost populations (GENELAND group D) contained only haplotypes (H33--H40) separated by at least 14 substitutions from the other haplotypes (H1--H32), but we did not detect this high degree of differentiation using microsatellites. A similar pattern of east-west differentiation was observed in the stonefly *Dinocras cephalotes* in the same study area \[[@CR62]\]. Specifically, two highly divergent haplotype groups were found for the *CO1* gene. However, contrary to our findings, Elbrecht *et al*. \[[@CR62]\] found haplotypes of these groups to be shared across populations of these groups and also detected ongoing gene flow using nuclear genes. As the stonefly has a terrestrial and more mobile life stage, these contrasting patterns may be explained by differences in mobility. The observed divergence between the two *G. fossarum* clade 11 groups (east-west) is likely the result of independent historic isolation in eastern and western refugia, as suggested by Elbrecht et al. (2014) for the stonefly species. Even though *G. fossarum* populations belonging to the two divergent groups seem to be isolated, it is likely that they can still interbreed when in contact. This is indicated by the low divergence (between 2.51 and 3.95 %) in comparison to previous estimates of divergence between cryptic species of *G. fossarum* \[[@CR10]\], also visible in the phylogenetic tree, when comparing intra- and inter-clade divergence*.* Further, Lagrue et al. \[[@CR63]\] found reproductive isolation between *G. fossarum* clades only when *CO1* sequence divergence was greater than 4 %.

We did, however, not only find strong differentiation between populations at a regional scale, but contrary to our hypothesis also at a local scale. Actually, most of the populations were highly isolated from each other (80 % of pairwise *CO1 F*~ST~ values were significant), sometimes even within 2 km in the same stream. The strong isolation is further supported by the pattern that many populations contained private haplotypes, which mostly differed by only a single mutation from the main haplotype of the respective population. This indicates long-term isolation because independent mutations were able to accumulate in the populations. In contrast to the overall strong isolation, we found little or no differentiation between some of the populations separated by over 40 km straight-line distances and even higher river distances (especially within GENELAND groups A and D) for the *CO1* gene. However, when analyzing a subset of these populations based on the more rapidly evolving microsatellites, we detected strong differentiation between all populations. Therefore, the subtle differentiation found between populations of group A is most likely not caused by ongoing gene flow; instead, low genetic variation is likely due to bottlenecks and the effects of genetic drift \[[@CR6], [@CR64]\].

In their study, Westram et al. \[[@CR24]\] concluded that differences in population structure between clade 11 and 12 could hint at interspecific differences in dispersal ability, life history or population size. The authors discussed that differences in differentiation levels reflect more likely species specific rather than being driven by geographic effects. However, the strong local differentiation we found in clade 11 resembled much more the pattern described for clade 12 than that of clade 11 \[[@CR24], [@CR36]\]. The discrepancies observed between the data of Westram et al. \[[@CR24]\] and our study suggest that geographic effects can strongly impact on the differentiation of populations and with that probably on the realized dispersal of populations within a species.

Based on the SHM \[[@CR37]\] our second hypothesis was that the populations are structured according to catchment boundaries and show an IBD pattern within catchments. The structuring of populations according to catchments has been shown for other crustaceans \[[@CR65], [@CR66]\] and an IBD pattern has been detected in *G. fossarum* \[[@CR12], [@CR34], [@CR36]\], specifically in clade 11 \[[@CR24]\]. Studying variance partitioning and using the sub-catchment boundaries in an AMOVA, we detected significant partitioning of genetic variance, but most of the variance (59.8 %) was found between populations within groups, indicating that populations were not primarily structured by sub-catchments, but were locally isolated within sub-catchments. We found a similar pattern when analyzing the correlation between genetic and geographic distances, as the correlation was significant, but very weak for both geographic distance metrics (straight-line and waterway distance), especially for waterway distances. The weak signature of sub-catchment-based structuring and the stronger correlation between genetic and straight-line distances in comparison to waterway distances indicate that processes other than low dispersal within streams are relevant to the genetic structure within this species. This inference is also supported by the GENELAND analysis, where groups A--E contained populations from at least two sub-catchments each and two populations in group B were geographically distant to all other populations of this group. Evidence for overland dispersal in aquatic invertebrates, including *G. fossarum*, was reported previously (e.g. \[[@CR36], [@CR67], [@CR68]\]), and it is assumed to occur by transport via vectors like birds \[[@CR69]--[@CR71]\], large mammals \[[@CR72], [@CR73]\], or humans \[[@CR13]\]. As some of the springs of the different sub-catchments in our study area are only separated by a few hundred meters and not by mountain ranges, overland dispersal is likely possible. However, as pairwise differentiation was high in general, overland dispersal seems to be more important for rare colonization events on evolutionary time scales rather than for recurrent dispersal at ecological time scales within generations. However, we did not expect the low correlation between waterway and genetic distances. A lack of IBD in other aquatic invertebrates is often attributed to particularly weak (e.g. \[[@CR6], [@CR74]--[@CR76]\]) or strong dispersal abilities (e.g. \[[@CR27], [@CR76]--[@CR78]\]). These results may also be explained by the presence of strong dispersal barriers between populations \[[@CR79]--[@CR81]\]. The slight IBD pattern we found in our study was mainly caused by the strong differentiation between the western- and easternmost populations, as they were separated by the greatest geographic distances and did not share any haplotypes. Otherwise, for many pairwise comparisons, we detected high genetic differentiation at small geographic distances and no differentiation over long distances. As we found strong population differentiation between populations, indicated by high and significant *F*~*ST*~ values, the lack of IBD cannot be caused by too strong realized dispersal. However, previous studies found IBD in different areas for *G. fossarum* (e.g. \[[@CR24]\]); accordingly, the dispersal ability should generally be sufficient to generate this pattern even though the realized dispersal appears to be very limited here. Therefore, it seems likely that other barriers to gene flow exist. These barriers could either be direct, like dams, weirs, or water reservoirs, or indirect if the conditions in connecting areas are unfavorable due to, for example, anthropogenic land use, organic pollution, acidification, large connecting rivers, or strong competition (e.g. \[[@CR20], [@CR27], [@CR36], [@CR68], [@CR74], [@CR76], [@CR82], [@CR83]\]). To analyze the influence of direct in-stream barriers, we used a subset of populations to determine if the number of barriers was correlated with the genetic distance between populations. Based on Mantel tests conducted for barriers of \>0.5 m and \>1 m, we did not observe a correlation with genetic distance, indicating that population isolation did not simply reflect the number of barriers. However, the influence of these barriers on realized dispersal could not be determined using the current marker and sampling scheme. To infer the influence of in-stream barriers more directly, individual barriers should be tested using more rapidly evolving genetic markers. Some of the sampling sites were separated from each other by water reservoirs, and, with one exception, the separated populations were significantly differentiated from all other analyzed populations, indicating a strong barrier effect of these water reservoirs. One exception at a first glance was in the Möhne sub-catchment, where populations were dominated by or even consisted exclusively of H1. However, as discussed previously, the diversity for the *CO1* gene in this region was too low to detect barrier effects. Therefore, we conclude that in-stream barriers, especially water reservoirs, can have a substantial impact on dispersal, yet they do not explain the general population structure observed in this study. With the exception of the water reservoirs that separate the clearly differentiated groups F and G from all other populations there were no obvious weir- or dam-related boundaries separating the groups.

In view of the tested parameters that cannot explain population structure sufficiently the question arises, which other factors ultimately underlie the apparent population structure. To determine if the structure was more influenced by colonization history (represented by the geographic position of the population) or by environmental factors that may differ between population groups, we conducted a discriminant analysis for the four main groups including variables for both possible drivers. The best predictor for the assignment of populations to the four GENELAND groups was longitude, indicating that the location of each site from east to west determines the affiliation to the *CO1* clusters. This clearly suggests that the observed population structure primarily reflects colonization history. If colonization events occur rarely and are initiated by relatively few individuals, this can lead to strong founder effects; genetic drift in small populations can lead to a strong loss of genetic diversity \[[@CR64], [@CR84]\]. The maintenance of this structure over time is probably promoted by a high degree of population isolation, leading to small effective population sizes and thereby enhancing the effect of random genetic drift. As *G. fossarum* is sensitive to many anthropogenic stressors, like organic pollution, a lack of oxygen, acidification \[[@CR15]\], and high ammonium concentrations \[[@CR14]\], isolated populations likely underwent drastic population declines over time, resulting in the loss of genetic diversity in this anthropogenically impacted area. Even though we could not identify a single anthropogenic factor influencing the population structure, we found populations to be highly isolated on local scales, clearly indicating that there are barriers to gene flow, which were not detectable with the methods used here.

Conclusion {#Sec12}
==========

In this study we found a considerably higher differentiation between populations of *G. fossarum* clade 11 in the human-impacted Sauerland area than was expected based on previous genetic analyses. The strong isolation was supported by two independent molecular marker systems, indicating that the realized dispersal was low in the study area. Also contrary to our hypotheses we found only a slight isolation by distance (IBD) pattern and structuring of populations according to river catchments but rather a strong geographic pattern (east-west differentiation). In view of published data, the dispersal ability of *G. fossarum* clade 11 specimens should be sufficient enough to create an IBD pattern. In the absence of this we conclude that there are barriers preventing gene flow partly, even between neighboring populations. Despite the likely effect of larger reservoirs on connectivity, we could not determine specific anthropogenic factors that directly influence the population structure. In fact, it was best predicted by the independent action of genetic drift at local sites after initial colonization. These effects are likely enhanced by the multitude of anthropogenic stressors, because *G. fossarum* is sensitive to many anthropogenic stressors. As *G. fossarum* clade 11 is widely distributed and is frequently found in the Sauerland area we conclude that its colonization ability over long time scales is good. The same holds likely true for its ability to establish new populations based on few colonizers, since long distance dispersal by animal vectors is unlikely to occur frequently \[[@CR3]\]. Even though we could not determine specific anthropogenic factors hindering gene flow, it is likely that the dispersal is influenced by these factors because of the exceptionally strong differentiation we found here in contrast to expectations based the findings of an earlier study (Westram et al. \[[@CR24]\] and the much broader distribution of clade 11 in comparison with the other *G. fossarum* clades. These findings highlight the importance to take regional factors into account when predicting the dispersal ability of species. Further more research is needed to determine the most important barriers restricting gene flow between populations to prevent further losses of genetic diversity and maintain an intact ecosystem.

Abbreviations {#Sec13}
=============

AMOVA, Analysis of molecular variance; *CO1*, cytochrome c oxidase 1 gene; FDR, False discovery rate; *F*~*ST*~, Fixation index; HWE, Hardy-Weinberg equilibrium; IBD, isolation by distance; SHM, Stream Hierarchy Model; VIF, variance inflation factor

Additional files {#Sec14}
================

Additional file 1:Ecological parameters for individual sampling sites. Group indicates the GENELAND group association of the sampling sites. Site is the abbreviation for each sampling site, and geographic coordinates are given in the Gauss--Krueger coordinate system. NA: no data were available for this site. (PDF 66 kb)Additional file 2:DNA salt-extraction protocol, modified from Sunnucks & Hales \[[@CR44]\]. (PDF 51 kb)Additional file 3:*CO1* haplotype information for each sampling site. Group indicates the GENELAND group association of the sampling sites, site is the abbreviation for the sampling site, and n refers to number of analyzed specimens per site. \# H gives the total number of haplotypes at a site and private H is the number of private haplotypes. H diversity is the haplotype diversity and H1 to H40 represent the different haplotypes. (PDF 67 kb)Additional file 4:**a:** Sampling locations of specimens used in the phylogeny. Symbols in the map correspond with symbols in **b** to indicate sampling locations of used sequences; the red square indicates the Sauerland area in which the specimens belonging to the main haplotypes were sampled. **b:** Neighbor-joining tree of main haplotypes and additional sequences belonging to clades 3, 10, 11, 12 and 13 (sensu Weiss et al. \[[@CR10]\]). The tree is drawn to scale, with branch lengths (next to the branches) in the same units as those of the evolutionary distances used to infer the phylogenetic tree (K2P method). (PDF 2114 kb)Additional file 5:Microsatellite loci information for the different sampling sites. Site is the abbreviation for the sampling site and n refers to number of analyzed specimens per site. AR is the rarefied average allelic richness over loci, and \# alleles is the number of alleles found for each locus. H~O~ is the observed and H~E~ the expected heterozygosity, and bold values indicate significant deviations from HWE. (PDF 44 kb)Additional file 6:Pairwise *CO1 F* ~ST~ values between all sampling sites. Sampling sites are indicated by the site abbreviations, and letters from A--G indicate GENELAND groups. Colors indicate within-group *F* ~ST~ values. Asterisks at site names indicate populations where microsatellite loci were also analyzed. Red and bold *F* ~ST~ values indicate significant values. Significance levels were adjusted according to the FDR. (PDF 67 kb)Additional file 7:Pairwise microsatellite *F* ~ST~ values between all sampling sites. Sampling sites are indicated by the site abbreviations, and letters from A--G indicate GENELAND groups. Colors represent within-group *F* ~ST~ values. Red and bold *F* ~ST~ values indicate significant values. Significance levels were adjusted according to the FDR. Below the diagonal, uncorrected *F* ~ST~ values are shown; above the diagonal, ENA-corrected *F* ~ST~ values are shown. (PDF 42 kb)Additional file 8:Percentage of populations predicted to belong to different GENELAND groups according to the discriminant analyses. (PDF 35 kb)

We thank Hannah Weigand, Jan Macher, Vasco Elbrecht, Markus Patschke, Vivienne Dobrzinski, Lisa Pöttker, Janis Neumann and Tobias Traub for assistance with sampling, Maria Gies and Martin Sondermann for providing some of the samples and Hannah Weigand, Vivienne Dobrzinski, Stefan Gerzmann, Markus Liebert, Anna Eckert and Julia van der Mond for help in the lab. We also gratefully acknowledge the help of the Ruhr-Verband and LANUV in providing the used ecological data. We further thank Ralph Tollrian, Hannah Weigand, Alexander Weigand, Jan Macher, Arne Beermann and Christian K. Feld for helpful discussions and three anonymous Reviewers from Axios Review for valuable comments. We thank the Bioedit copy editing service for correcting language in an earlier version of the manuscript.

Funding {#FPar1}
=======

This work was funded by a junior research group grant of the Kurt Eberhard Bode Foundation to FL.

Availability of data and material {#FPar2}
=================================

Sequences of the 40 haplotypes are available in GenBank under the accession numbers KX065366 -- KX065405. Individual Sequences and microsatellite data can be made available upon request.

Authors' contributions {#FPar3}
======================

MW and FL conceived and designed the study. MW collected the samples, performed the laboratory work and the data analysis. MW wrote the manuscript together with FL and both authors approved the final version of it.

Competing interests {#FPar4}
===================

The authors declare that they have no competing interests.

Consent for publication {#FPar5}
=======================

Not applicable.

Ethics approval and consent to participate {#FPar6}
==========================================

No special permissions were required for sampling *G. fossarum*. We obtained special permissions to collect specimens at sites in nature reserves.
