Individual small in‐stream barriers contribute little to strong local population genetic structure five strictly aquatic macroinvertebrate taxa

Abstract Water flow in river networks is frequently regulated by man‐made in‐stream barriers. These obstacles can hinder dispersal of aquatic organisms and isolate populations leading to the loss of genetic diversity. Although millions of small in‐stream barriers exist worldwide, their impact on dispersal of macroinvertebrates remains unclear. Therefore, we, therefore, assessed the effects of such barriers on the population structure and effective dispersal of five macroinvertebrate species with strictly aquatic life cycles: the amphipod crustacean Gammarus fossarum (clade 11), three snail species of the Ancylus fluviatilis species complex and the flatworm Dugesia gonocephala. We studied populations at nine weirs and eight culverts (3 pipes, 5 tunnels), built 33–109 years ago, mainly in the heavily fragmented catchment of the river Ruhr (Sauerland, Germany). To assess fragmentation and barrier effects, we generated genome‐wide SNP data using ddRAD sequencing and evaluated clustering, differentiation between populations up‐ and downstream of each barrier and effective migration rates among sites and across barriers. Additionally, we applied population genomic simulations to assess expected differentiation patterns under different gene flow scenarios. Our data show that populations of all species are highly isolated at regional and local scales within few kilometers. While the regional population structure likely results from historical processes, the strong local differentiation suggests that contemporary dispersal barriers exist. However, we identified significant barrier effects only for pipes (for A. fluviatilis II and III) and few larger weirs (>1.3 m; for D. gonocephala). Therefore, our data suggest that most small in‐stream barriers can probably be overcome by all studied taxa frequently enough to prevent fragmentation. However, it remains to be tested if the strong local differentiation is a result of a cumulative effect of small barriers, or if larger in‐stream barriers, land use, chemical pollution, urbanization, or a combination of these factors impede gene flow.


Introduction 45
Rivers in the Anthropocene are in crisis due to various human activities (Jackson et al., 2001;46 Vörösmarty et al., 2010). In particular, natural water flow has been altered substantially over the past 47 ddRAD data analysis 153 For a better overview, all important analysis steps are summarised in the workflow in Figure 3.  processing of ddRAD data for each sequencing library was conducted as described in Weiss et al. (2018). 155 After pre-processing, denovo_map.pl of Stacks v.1.34 (Catchen et al., 2013) was used to identify loci 156 and genotypes. Stacks was run with eight different settings per species to identify the optimal parameter 157 settings (Paris, Stevens, & Catchen, 2017, see Appendix S1 for details). Stacks results were exported 158 from stacks databases with export_sql.pl (minimum stack depth 8). For D. gonocephala, aneuploid 159 cytotypes are known besides diploids, with chromosome numbers similar to triploids (de Vries, 1986). 160 Therefore, individual ploidy levels based on the expected allelic coverage was estimated using the R-161 script ploidyCounter.R (Rozenberg, https://github.com/evoeco/radtools/, see Weigand et al. 2018 for 162 details). However, this ploidy estimation is indirect and aneuploids with high chromosome numbers 163 cannot be distinguished from real triploids. Therefore, specimens with coverage distribution plots as 164 expected for triploids are called 'potential triploids' hereafter. Direct estimation of ploidy levels, e.g. 165 via flow cytometry was not possible from the ethanol fixed material. All analyses were also conducted 166 excluding potential triploids, but as they produced similar results only analyses including all specimens 167 are presented here. 168 Further analyses were performed using Snakemake workflows (Köster & Rahmann, 2012), 169 combining stacks2fasta.pl (Macher et al., 2015) and several in-house python and R-Scripts for data 170 reformatting, filtering and population genetic analyses. Parameter settings were similar for each species, 171 except for Ancylus species, where parameters were adjusted when possible to account for polyploid 172

genomes. 173
First, Stacks parameter and locus filtering settings were optimised, as described in detail in 174 Appendix S1. General population structure was analysed with individuals divided into populations 175 according to barrier sites, i.e. pooling specimens from individual sampling sites at a barrier site to one 176 population ( Figure 3, step 6). The following locus filtering settings were chosen: 1-12 SNPs/locus (only 177 one used), minor allele frequency 1%, locus present in at least 80% of individuals per population and 178 95% of all individuals. Further, specimens with too few reads, or more than 20% missing data were 179 excluded. Basic population genetic statistics, e.g. observed heterozygosity (HO), observed gene diversity 180 (HS), overall gene diversity (HT) and overall FST and FIS were calculated using the R-package hierfstat 181 (Goudet, 2005)  and k = 1-17 for G. fossarum and D. gonocephala, and analyses were run with 30 replicates and 100,000 188 iterations per replicate. In addition, Neighbour-Net networks (Bryant & Moulton, 2004) were calculated 189 using SplitsTree v. 4.14.5 (Huson & Bryant, 2006). Pairwise FST values (after Weir & Cockerham, 1984) 190 between barrier sites were calculated and significance was tested by bootstrapping over loci (10,000 191 replicates; 0.025/0.975 confidence intervals) using the R-package hierfstat. HO and allelic richness (AR) 192 in each population were calculated using the R-package diveRsity (Keenan et al., 2013). Further, the 193 divMigrate function (Sundqvist et al., 2016) of this R-package was used to assess directional relative 194 migration rates and to detect asymmetries in gene flow using Nei's GST (Nei, 1973). To evaluate 195 correlation of genetic and geographic distances (isolation by distance, IBD), Mantel tests were 196 conducted using the R-package vegan (Oksanen et al., 2019). As a measure of genetic distance, pairwise 197 FST values were calculated between all individual sampling sites (S). For geographic distances, either 198 straight-line or waterway distances were used. Both distances were calculated using QGIS v. 2.14.14 199 (http://qgis.org) with the same stream map used for the visualisation of sampling sites and population among adjacent sites, loci were exported separately for each barrier site with the same filter settings 205 used for previous analyses, except the minor allele frequency was set to 5% and loci had to be present 206 in 90% of specimens to account for the smaller sample sizes. For barrier analyses, two population 207 definitions were used. First, all analyses were conducted with barrier sites as individual populations, i.e. 208 specimens from individual sampling sites up-and downstream of a barrier were pooled (S upstream/S 209 downstream). Second, all individual sampling sites (S) were treated as different populations (see Figure  210 3). Basic population genetic statistics, sNMF (k = 1-6), HO and AR per population and pairwise FST were 211 calculated as before. The number of private alleles in each population was calculated using an in-house 212 python script. To account for differences in sample size, larger populations were randomly subsampled 213 30 times to match the smallest population size and the mean number of private alleles among replicates 214 was calculated. DivMigrate was used to calculate migration rates between sampling sites and to assess 215 asymmetric migration patterns. A large number of loci (>1000) can partly compensate for the low 216 sample sizes (Willing et al., 2012), but results for populations with n < 5 should be interpreted with 217 caution (Sundqvist et al., 2016). To test the hypothesis that small in-steam barriers generally increase 218 fragmentation compared to reference sites and that weirs have a stronger effect than pipes or tunnels 219 (both culverts), linear mixed-effect models (LMMs) were run using the R-packages lme4 (Bates et al.,220 2015) and lmerTest (Kuznetsova et al., 2017). First, mean differentiation (FST) for reference sites and 221 across barriers were calculated per barrier and species. Different LMMs were run for each species with 222 differentiation (FST) as dependent variable, presence of barriers between compared sites as fixed effect 223 and barrier site as random effect. In a second LMM analysis, barriers were subdivided into barrier types, 224 i.e. weirs, pipes, and tunnels (and no barrier) and compared per species. 225 226

Population genetic simulations 227
To simulate expected population differentiation caused by migration barriers over time, Nemo 2.3.51 228 (Guillaume & Rougemont, 2006) was used. As there is only little information available about the 229 parameters needed to be given in the simulation, various barrier strengths and biological species traits 230 were applied, including different mating systems (random mating for G. fossarum and hermaphrodite 231 for A. fluviatilis and D. gonocephala), population sizes, extinction rates, selfing rates for hermaphrodites 232 and dispersal and barrier models (see Figure 3 and Appendix S2 for details on parameter space). The 233 dispersal and barrier models were designed to reflect stream dispersal (linear, influenced by flow 234 direction) and our sampling scheme (four equidistant sampling sites). To prevent edge effects (no further 235 dispersal possible for edge populations in the simulation), we simulated six populations but only 236 analysed the four target populations. Simulations were run for 250 generations. The first 100 generations 237 were simulated with different dispersal models but without a barrier to separate dispersal model and 238 barrier effects. After 100 generations, a barrier was introduced for the next 150 generations using three 239 models (complete barrier, upstream barrier with unchanged, or reduced downstream dispersal, see 240 All five macroinvertebrate species showed strong regional and local population structure as detected by 256 all methods. 257

Population structure of Gammarus fossarum 258
We found G. fossarum in nine streams at 11 barrier sites (Table S1, Tables S3 and S4 for summary  259 statistics for all species). PCA ( Figure S2) and sNMF analyses indicated strong population structure. In 260 the sNMF analysis, cross-entropy (CE) was lowest for k = 9 ( Figure S3A). The nine clusters 261 corresponded to the nine streams, also visible in the Neighbour-Net analysis ( Figure 4A and B). Further, 262 the network showed superordinate clustering into three phylogenetically distinct groups in concordance 263 with sampling site topography. Pairwise FST values were generally high (mean 0.47), especially between 264 network-groups and all values were significant, including values for the comparison between barrier 265 sites in one stream, ranging between 0.01 and 0.74 ( Figure 4C). Likewise, estimated migration rates 266 were very low in most cases (mean 0.05), with the exception of within-stream comparisons ( Figure 4D, 267 Table S5A). However, for the two within-stream comparisons, migration rates between VR6 and VR17 268 were much lower in both directions than those between QB17 and QB23 and were significantly higher 269 in the downstream than in the upstream direction. In 71% of the comparisons, significant asymmetric 270 gene flow was detected. were potentially triploid. PCA ( Figure S2) and sNMF cluster analysis both indicated strong population 275 sub-structuring. CE was lowest for K = 13 ( Figure S3B). The population structure could be explained 276 best by the stream origin and ploidy level ( Figure 4E). At sites with both diploids and potential triploids 277 (QB24, QB17 and QB23 and VR9), individuals were divided into two clusters or showed high 278 intermixture between clusters. The Neighbour-Net showed a star-like pattern containing 11 clusters, 279 each representing a single stream ( Figure 4F), but sub-clusters within streams according to ploidy level 280 were visible. FST values were lower than those for G. fossarum (mean: 0.31, range 0.08 to 0.51) but 281 significant for all pairwise comparisons including within stream comparisons ( Figure 4G). In general, 282 estimated migration rates were higher than those for G. fossarum (mean: 0.16) and migration rates for 283 most comparisons were significantly asymmetric (73%; Figure 4H, Table S5B). Similar to G. fossarum, 284 migration rates were highest between sites within a stream and migration was significantly higher in the 285 downstream direction.  (Table S4). For A. fluviatilis I, PCA ( Figure S2) and sNMF indicated strong 291 genetic structure. CE values were low k = 7-10 ( Figure S3C), with the lowest mean CE over repetitions 292 for k = 8. The clusters mainly corresponded to streams, but more specimens showed higher admixture 293 to additional clusters ( Figure 5A). These results corresponded well with the network results ( Figure 5B), 294 where a similar group structure consistent with streams was visible, with the exception of specimens 295 from QB12 and QB20, which were split into two groups according to the barrier site. In general, pairwise 296 FST values were lower than those in G. fossarum and D. gonocephala (mean: 0.18, range: 0 -0.29, 297 Figure 5C), but all were significant with the exception of the within-stream comparison QB11 and 298 VR20. Migration rates were between the two previously described species (mean: 0.11), with 80% of 299 comparisons showing significant asymmetry ( Figure 5D, Table S5C). 300 For A. fluviatilis II, PCA ( Figure S2) and sNMF both indicated strong population structure 301 within the species. The most probable number of clusters in the sNMF analysis was k = 8 ( Figure S3D). 302 The clusters were less clear than those for the other species, as many specimens showed high amounts 303 of shared ancestry to different clusters. The main factor explaining the groups was stream origin ( Figure  304 5E), but in some cases, main cluster were shared among different streams or additional cluster were 305 found within one stream. A similar structure was observed in the Neighbour-Net analysis. However, 306 single specimens with high admixture in the sNMF plot, grouped together with specimens sampled in 307 other streams ( Figure 5F). Overall, pairwise FST values indicated relatively low but significant 308 differentiation between all sampling sites, except for QB23-QB24 (mean: 0.08, range: 0.003 -0.15, 309 Figure 5G). Estimated migration rates were higher than those for the other species (mean: 0.17) and less 310 asymmetry was detected (58% of comparisons, Table S5D). 311 A. fluviatilis III was mainly found in one stream with two barrier sites (VR6 and VR17). At 312 VR9, four individuals were identified as A. fluviatilis III/I hybrids (75%/25%) according to Weiss et al. 313 (2018) and these were included here. In the PCA ( Figure S2) and sNMF analysis, sub-structure was 314 detected, with four clusters best explaining the population structure ( Figure S3E). One cluster consisted 315 of VR9, the second of VR17 specimens and the other two of VR6 specimens, with shared proportions 316 for some individuals between clusters ( Figure 5I). A similar pattern was observed in the Neighbour-Net 317 analysis ( Figure 5J). Pairwise FST values were all significant but lower between VR6 and VR17 ( Figure  318 5K), where also higher migration rates were detected ( Figure 5L, Table S5E). 319 For all species, we detected a correlation between geographic and genetic distances, which was 320 strongest for A. fluviatilis II and weakest for D. gonocephala ( Figure S4). Genetic diversity (i.e., HO and 321 AR) differed among sites for all species, but with no consistent pattern over all species (Table S6). 322 323

Effects of in-stream barriers on population structure 324
Filtering loci separately for each barriers site, resulted in in a greater number of variable loci and mostly 325 increased genetic diversity per barrier site in comparison to the total dataset (Table S7). We could not 326 detect distinct barrier effects across species when comparing FST values between reference and barrier 327 separated sites. Mean differentiation was low and mostly insignificant regardless of the presence of 328 barriers or the barrier type (weir, pipe or tunnel). Consequently, none of the LMMs indicated a general 329 effect of barriers on differentiation (Appendix S3). When studying specific barriers in detail, we found 330 no evidence for a barrier effect on population structure for seven of the nine weirs and six of the eight 331 culverts in any species. This means that populations at barrier sites were either i) not differentiated at 332 all, ii) differentiation was similar among reference and barrier separated sites, iii) lower for some of the 333 comparisons of barrier separated sites compared to reference sites, or iv) that only individual sampling 334 sites were constantly differentiated to the other sites from the same barrier site. All these four 335 differentiation patterns were found for all species. Differentiation between individual sites was most 336 prominent for specimens sampled in affluents ("N"; Figure 6), despite the lack of obvious barriers. At 337 three weirs, reduced upstream migration rates were detected for A. fluviatilis I and II (QB11, QB12 and 338 QB17) ( Figure S5). However, we found no indication for a barrier effect in any of the other analyses. A 339 similar reduction was observed for VR6a for D. gonocephala but here the sample size was too small at 340 upstream sites ( Figure S5). 341 We only detected distinct barrier effects associated with weirs for D. gonocephala at QB24 and 342 QB11. In both cases, differentiation across barriers was larger than that between reference sites ( Figure  343 6). The effect was stronger for QB24 than for QB11, consistent with the sNMF results, where two 344 clusters were the most probable for QB24. Specimens sampled upstream of the barrier only belonged to 345 the first cluster, while specimens assigned to both clusters, as well as intermixed ones, were found 346 downstream of the barrier. Specimens showing high membership to the second cluster were all classified 347 as potential triploids. At both barriers, migration rates in up-and downstream directions were lower 348 across barriers than among reference sites. For QB24, all rates across the barrier were significantly 349 asymmetric, with stronger downstream dispersal. Also, genetic diversity was higher downstream of 350 barriers at both sites, indicated also by higher amounts of private alleles (Tables S8 and S9). 351 A significant barrier effect was detected for A. fluviatilis II at VR20 and for A. fluviatilis III at 352 VR6b. While VR6b was a pipe, VR20 was a combination of a pipe and tunnel with an additional steep 353 drop within and two streams were joined within the tunnel. At VR20 S4, we only found one specimen 354 belonging to A. fluviatilis II, which was excluded from analyses of FST and migration rate. The other 11 355 specimens at this site were identified as A. fluviatilis I, which was not detected upstream of the culvert. 356 All individual site comparisons across the culvert in A. fluviatilis II indicated significant differentiation, 357 including comparisons of both N1 and N2 with S2, which were separated by a pipe and part of the 358 tunnel, while differentiation among reference sites was not significant. This pattern was also detected in 359 the divMigrate analysis. Upstream migration across the culvert to both ends was highly reduced, while 360 downstream migration rates were similar or even higher than those among reference sites ( Figure S5). 361 In the sNMF analysis, indications for a barrier effect were found with the lowest CE for three clusters. 362 Specimens from both N1 and N2 were split into two clusters (k1: n = 10, k2: n = 6), S1 specimens were 363 assigned to k2 or the third cluster (k2 = 6, k3=5) and specimens at S3 were assigned to all three clusters, 364 indicating that downstream dispersal was possible. AR was lowest at S2; otherwise, genetic diversity 365 was similar (Tables S8 and S9). 366 At barrier site VR6, the longer pipe of this barrier site (VR6b) caused population differentiation 367 in A. fluviatilis III, as all comparisons across the pipe were significant and FST values were higher than 368 those between reference sites ( Figure 6). These differentiation patterns were supported by the sNMF 369 analysis, where three clusters were detected. One cluster was only found in S3&4, while the other two 370 were mainly found in S1&2 and less frequently in S3&4, where assignment to multiple clusters was 371 observed. While migration rates between sites were generally reduced, upstream migration across the 372 second pipe (VR6b) was even stronger reduced to nearly zero. By contrast, migration rates across the 373 first pipe (VR6a) were even higher than those among reference sites ( Figure S5). HO was similar at all 374 sampling sites, while AR and the number of private alleles were higher at S3&4 than at other sites 375 (Tables S8 and S9). The simulation of a complete barrier revealed the quick emergence of strong population differentiation, 381 as determined by increases in FST values and decreases in migration rates over time ( Figure S6, Table  382 S10). In general, barrier effects were similar, independent of the distance between populations. In 383 smaller populations, fluctuations over time were larger, leading to low but significant population 384 differentiation between reference populations not separated by a barrier. However, after only a few 385 generations, FST values for populations separated by a barrier were consistently higher than those 386 between reference populations. For larger populations (i.e. patch capacity 2000), it took up to 50 387 generations to detect differentiation and FST values remained low but followed the same general pattern 388 ( Figure S6). Here, only FST values between populations separated by a barrier were significant. The 389 higher the extinction rate, the higher the FST values between barrier-separated populations and also in 390 larger populations, differentiation was detected earlier. For smaller populations with higher extinction 391 rates, significant differentiation was detected without a barrier but remained lower than for barrier-392 separated populations. These patterns were visible for all tested migration matrices, with asymmetric 393 matrices leading, on average, to slightly higher FST values. However, while a reduction of migration 394 rates over time was detected with the divMigrate approach, the simulated asymmetry was not reflected 395 in consistent differences between up-and downstream migration rates. Applying different selfing rates 396 to simulate hermaphrodites, FST values increased with increasing selfing rates, but patterns were 397 generally the same as those for populations with obligate sexual reproduction. Even though complete 398 isolation between populations on both sides of the barrier was not directly detected, a reduction in 399 migration rates after barrier introduction was quickly (<10 generations) detected for all migration 400 scenarios and decreased further with time. The reduction was more easily detected for larger sample 401 sizes, higher extinction rates and smaller populations. 402 In general, barriers impeding upstream but not downstream dispersal did not affect the population 403 structure and did not result in reduced inferred migration rates ( Figure S6, Table S10). At some time 404 points for small populations (particularly in combination with higher extinction and/or selfing rates), 405 FST values were slightly higher and significant for barrier-separated populations than for reference 406 populations, but this pattern fluctuated over generations and did not increase with time (similar for 407 migration rate estimates). Reducing downstream dispersal rates led to an intermediate pattern between 408 those described previously. Not all parameter combinations were simulated, but downstream dispersal 409 rates had to be reduced substantially to 1 %, or depending on the other parameters 2 %, to influence 410 population differentiation. The effects were stronger (and sometimes only detectable) for smaller 411 population sizes, higher extinction or selfing rates, and were strongest for a combination of these 412 parameters (i.e., if population sizes were small and extinction rates were high). If the selfing rate was 413 increased from 0.3 to 0.6, the effect on the population structure was generally high but barrier effects 414 were less detectable because fluctuations between populations and over generations were higher. We 415 obtained similar results for migration rates, but correct asymmetry was not reliably detected. HO and AR 416 were, in general, not affected by the barrier model ( Figure S7). Both AR and HO decreased stronger over 417 time for smaller population sizes and higher extinction proportions for all dispersal models (with no or 418 slow decreases for larger population sizes, i.e., 1000 and 2000). However, the decrease was constant 419 over time, independent of the introduction of a barrier into the system and similar in all populations.

General effects of small in-stream barriers on macroinvertebrate taxa 424
For all five stream species, our population genomic data revealed strong population differentiation at 425 local scales. However, in contrast to our hypothesis, the tested individual small barriers did only increase 426 genetic differentiation compared to control sites in few cases. No correlation between barrier strength 427 (weir, pipe or tunnel) and differentiation was found. Therefore, we found little evidence that the tested The simulation results consistently showed that barriers preventing up-as well as downstream 438 dispersal lead to strong population structure relatively quickly within the age range of the barriers tested. 439 Other studies simulating the effect of barriers on population structure found that barrier effects will only 440 be detectable with divMigrate, complete barriers were detectable as well by a general reduction in migration rates between 450 barrier-separated sites. However, in general, migration rates were overestimated, suggesting that 451 inferred low migration rates in real data probably indicate real barrier effects. 452 To conclude, simulations showed that even minor migration (2 % or 5 % depending on population 453 size, extinction and selfing rate) in one direction can be sufficient to counteract barrier-related 454 population differentiation. If populations are large (patch capacity >2000) and no extinction rates are 455 simulated, even 1 % downstream migration can be enough to counteract differentiation in random 456 mating populations. Accordingly, reduced migration due to barriers is expected to be detectable if 457 migration rates are low, but detection will be difficult or impossible for very young (<10 years) barriers, 458 for barriers only reducing upstream dispersal and for large effective population sizes. While no information is available on the probability or frequency of zoochorous dispersal for the taxa 474 studied here and in general more information is needed, e.g. on the quantitative contribution of zoochory 475 to dispersal (Coughlan et al., 2017), it should be considered as a possible mechanism counteracting 476 barrier effects. 477 With respect to the individual species, the barrier effects were weakest for G. fossarum clade 11, 478 as expected. This indicates that dispersal ability is high enough to overcome the small in-stream barriers 479 analysed here (see also , congruent with findings for G. fossarum clade 12 (Alp 480 et al., 2012). It is also possible that the tested barriers were not old enough to create detectable isolation 481 patterns (Monaghan et al., 2001), that the reduction in migration was not strong enough to impact Our prediction that barrier effects would be moderate for D. gonocephala and strongest for A. 489 fluviatilis was not supported by the data. Only two weirs and two culverts had species-specific barrier 490 effects. The weirs influencing D. gonocephala were classified as the second (QB24) and third (QB11) 491 most severe according to height, steepness and smoothness. They were both >80 years old. For other 492 barriers of the same age, no population structure was detected as well as for QB27, which was even 493 higher than QB24 with a smoother ramp slope. The difference between barriers of the same age might 494 be explained by differences in population size or by important barrier characteristics that had not been 495 measured. The detected barrier effect was weaker for QB11 than for QB24 which was not as high and 496 had a moderate slope with drops at the start and the end. The two drops might have strengthened the 497 effect in comparison to other barriers with similar slope. In general, our results indicate that D. 498 gonocephala maintains gene flow across most small weirs, but we also found evidence that weirs of the 499 size and shape tested can impede migration. This holds true especially for those classified as more or 500 most severe. Thus, larger barriers, especially larger drops, will probably impact on dispersal for this 501 flatworm species. Contrary to our expectations, A. fluviatilis migration was not stronger affected by 502 weirs. At three weirs, upstream migration rates across barriers were lower than those between reference 503 sites, but no effect was detectable with any other method. To finally conclude that small weirs do not 504 influence migration in A. fluviatilis, it would be important to increase the sample size and to evaluate 505 the influence of polyploidy on analyses. In addition, especially larger drops should be investigated as 506 they will probably pose a greater barrier to dispersal than steep ramps, if they cannot be crossed by 507 crawling or passive zoochorous dispersal. 508 For culverts, patterns were more consistent with our expectations. We only detected effects for A. 509 fluviatilis II and III at two culverts predicted to have the strongest impact (VR20 and VR6b). At VR20, 510 only A. fluviatilis II was found upstream of the culvert, while A. fluviatilis I was additionally found at 511 downstream sites. While downstream dispersal seemed to be possible through this culvert, upstream 512 dispersal was highly limited. A similar pattern was found for A. fluviatilis III at VR6b, even though here, 513 inferred migration in both directions was reduced also between reference sites, but stronger across the 514 pipe. For the other species, we could not attribute observed population structure between sampling sites 515 to culverts, or the sample size upstream of the barrier was too small (D. gonocephala). These results 516 indicate that all studied species can probably disperse effectively through tunnels of up to 120 m length, 517 but pipes < 25 m can act as strong dispersal barriers probably due to their smooth internal structures and 518 high flow velocity (David et al., 2014), or drops at the end, preventing upstream movement (Vaughan,519 2002). Also for culverts it is possible that barrier effects went undetected due to, e.g. large Ne. But in 520 general our results were in accordance with the expectations that tunnels present less severe barriers 521 than pipes, and it seems possible that the studied taxa can overcome them actively. However, with 522 respect to pipes our results suggest that these can pose severe dispersal barriers even when relatively 523 short. Therefore, it would be important to focus on this kind of culverts in further studies. 524 In general, single small barriers did not have a large impact on population structure, yet it remains 525 to be determined whether cumulative effects exist and population structure increases when many small 526 barriers occur within short distances. Future work should focus on testing and adapting novel and 527 promising indices to identify population fragmentation in streams developed based on microsatellites 528 (Prunier et al., 2020) to SNP-based data including also possible asymmetric gene flow and focus on 529 barriers at the upper level of the severity gradient tested here, for weirs as well as for culverts. 530 531

Regional subdivision 532
While effects of individual in-stream barriers was minor, genetic differentiation between streams was 533 high for all species despite the small spatial scale of the study. Based on life history traits, we 534 hypothesised that population structure to be weakest for In D. gonocephala, we detected fairly high population differentiation, but in contrast to G. fossarum, no 563 strong geographic structure was found, supported by the low IBD pattern. Populations separated by only 564 approximately 2 km in the same stream already showed significant population differentiation. 565 Differentiation was also detected for 11 to 18 km but was not correlated with waterway distance within 566 this range, suggesting that additional factors other than waterway distance influence effective migration 567 in this species. Local adaptation could shape the population structure, as reported for a population in the 568 same study area in response to high copper concentrations , in which also high 569 population differentiation was found. Apart from this, little is known about effective dispersal in D. 570 gonocephala, but planarians are generally been regarded as weak dispersers (e.g., Rader et al., 2017). 571 The overall lower FST values in comparison to G. fossarum could indicate that D. gonocephala is a better 572 disperser. However, it is more likely that differences in FST reflect the greater phylogenetic signal in G. 573 fossarum, represented by three deep phylogenetic lineages as described above. In D. gonocephala, 574 particularly strong population structure in some cases was associated with the presence of potentially 575 triploid individuals. Potential triploids originating from different streams did not cluster together but 576 were closely related to specimens in the same stream, indicating that polyploidy evolved independently 577 several times. However, results concerning possible polyploid individuals have to be interpreted with 578 caution, i) because it is difficult to assess how analyses are affected for example by allelic dosage supported by the sNMF analyses. Additionally, most within-stream comparisons showed low but 600 significant differentiation, even within a few hundred meters. In general, our findings suggest that 601 effective dispersal is low between streams and differentiation within streams can occur over short 602 distances (see also Macher et al., 2016). In some cases, we detected genetic structure even within barrier 603 sites, while in other cases, levels of shared ancestry between clusters were higher than those for G. 604 fossarum or D. gonocephala. This may reflect higher occasional gene flow between streams, a younger 605 recolonisation history, or could have been influenced by polyploidy and the still high heterozygosity in 606 the data set. For a better understanding of the determinants of population structure, it is important to 607 better characterise homoeologous loci and to determine ploidy levels of individuals. However, the 608 consistent overall patterns across species suggest that our results are reliable and population isolation 609 despite differences in dispersal capacity indicates that strong dispersal limitations exist in the area, which 610 are not identified, yet. 611 612

Conclusions 613
We detected strong genetic isolation among populations of five hololimnic species between streams. 614 While for some individual barriers and species an effect was found, our genomic ddRAD data suggest 615 that single small weirs and culverts are probably not the cause for the detected strong fragmentation at 616 few kilometres distances. It remains to be tested, if cumulative effects of small barriers could have 617 caused the disruption of gene flow between populations, or if other factors, such as larger in-steam 618 barriers, land use, chemical pollution, urbanisation, or a combination of these, led to the observed 619 population structure. In general, our combination of genomic markers, population genetic simulations 620 and a controlled sampling design with distinct reference populations is a suitable tool to infer major 621 drivers of regional and local population structure. ddRAD data for G. fossarum and D. gonocephala will be available at NCBI BioProject upon publication. 634 ddRAD data for A. fluviatilis I, II and III will be available at NCBI BioProject with accession number 635

Conflict of interest 638
The authors declare that they have no competing interests.  Table 1. Barrier characteristics of culverts (VR) and weirs (QB), sorted by expected severity from high to low. 867 First three VRs were categorized as pipes, all other were tunnels. 868   colours. Barriers are ordered according to severity, starting with the strongest weir and strongest 902 culvert. Combi: differentiation across both barriers at one barrier site. Barriers where an effect was 903 detected are indicated with an asterisk. 904