The tale of springs and streams: how different aquatic ecosystems impacted the mtDNA population structure of two riffle beetles in the Western Carpathians

Introduction

Prolonged isolation of populations influences their genetic diversity and can be considered as the main force shaping genetic structure of aquatic species in Europe (Bálint et al., 2011; Alp et al., 2012; Theissinger et al., 2012). Isolated populations of aquatic invertebrates may be found within geomorphological units/subunits of many mountain areas (Engelhardt, Haase & Pauls, 2011; Davis et al., 2013; Mamos et al., 2016; Čiamporová-Zat’ovičová & Čiampor Jr, 2017; Šípošová, Čiamporová Zat’ovičová & Čiampor Jr, 2017; Copilaş-Ciocianu et al., 2018).

The W Carpathians are considered a biodiversity hotspot for a wide range of aquatic and terrestrial taxa (Neumann et al., 2005; Kotlík et al., 2006; Theissinger et al., 2012; Vörös et al., 2016; Copilaş-Ciocianu et al., 2017; Juřičková et al., 2017). However, the biodiversity of the W Carpathians is still underexplored, especially in terms of genetic diversity and population structure of aquatic species. In this context, studies upon the phylogeography of aquatic biota should be more focused on springs, or more generally, headwaters that are now heavily understudied compared to other aquatic biotopes.

Springs support unique macroinvertebrate communities that are found nowhere else in a catchment (Lewin et al., 2015). They are characterized by chemical, physical, and trophic constancy over several geological periods (Minshall & Winger, 1968; Odum, 1971; Butler & Hobbs, 1982; Cushing & Wolf, 1984; Glazier & Gooch, 1987; Pringle et al., 1988; Gooch & Glazier, 1991; Orendt, 2000; Wood et al., 2005; Meyer et al., 2007), which in turn provided a stable environment for aquatic invertebrates during adverse climatic conditions (Malicky, 2006; Ujvárosi et al., 2010). Springs function as ecotones between the surface and underground waters, which makes them an ecologically significant habitat (Gibert, 1991).

Accordingly, the main objective of our study is to compare genetic population structure and diversity patterns of two aquatic beetle species of Elmidae family, Elmis aenea (PWJ Müller, 1806) and Limnius perrisi (Dufour, 1843) in springs and streams of the W Carpathians. Said species are relatively closely related and commonly co-occur, yet in terms of population genetics represent a generally understudied family of freshwater beetles. Limited dispersal abilities, high habitat specificity, and more or less fragmented distribution make Elmidae an ideal taxon for studying genetic diversification through many geographic regions. Both studied species are rheophilic, oligo-stenotherm, and typical inhabitants of epirhithral streams at higher altitudes (Moog & Jäch, 1995; García-Criado, Fernández-Aláez & Fernández-Aláez, 1999). They are relatively widespread, which guarantees the detection of possible gene flow among geomorphological units/subunits of W Carpathians.

Our study aims to answer the following questions: (a) Are the spring subpopulations genetically more variable when compared to subpopulations in the streams? (b) Does genetic structuring of populations reflect population size change? (c) Are there interspecific differences in the population genetic structure among these related beetle species?

Materials & Methods

Study area

The Carpathians form an arc of mountains stretching across Central and Eastern Europe, with its main geomorphological units being the Western and Southeastern Carpathians (Kondracki, 1989). For this study, we focused on the W Carpathians representing the northernmost segment of the Alpine-Carpathian mountain chain (Fig. 1A). The W Carpathians reach medium altitudes (ranging from 500 to 1,300 m a.s.l.), only a few of their ranges exceed 1,500 m a.s.l.; geologically the mountain system is characteristic by interactions of rock folding and horizontal shifts (Bielik, 1999). Differences in altitude and in the geomorphological relief determine the precipitation in the area. In general, W Carpathian rivers have a rain-snow regime with floods in spring and summer.

Investigated localities (36 springs and 37 streams) are situated mainly on the territory of the Slovak Republic, partially in the Czech Republic and Poland; in the geomorphological units/subunits of the Inner and Outer Western Carpathians. The exceptions are Vihorlat Mts (VM) being part of the Inner Eastern Carpathians and Poloniny Mts (PM) belonging to the Outer Eastern Carpathians (Fig. 1B, Table S1).

Results

The distribution of E. aenea and L. perrisi suggests statistically significant different habitat preferences between these species in the W Carpathians (p < 0.0001; Fisher’s exact test). E. aenea has a rather wide distribution in karst springs (31 sites), while it is less widespread in streams (16 sites). On the contrary, L. perrisi was found only in eight springs, but in 30 streams. L. perrisi was also found in four streams of VM (Inner Eastern Carpathians) and in one stream of PM (Outer Eastern Carpathians), while E. aenea was not recorded in these geomorphological subunits. Both species co-occurred only in three springs and 14 streams from the total of 73 sites sampled in the W Carpathians (Table S1). There was also a significant effect of altitude on the presence of both species registered, but the dependence between altitude and habitat type (spring and stream) has not been demonstrated. The altitude range of both elmid species is shown by boxplots in Fig. 1C.

The haplotype distribution within investigated area shows that local subpopulations of the two elmid species are dominated, each, by the same widespread haplotype (Figs. 2A, 3A). Haplotype networks (Figs. 2B, 3B) also showed the similar haplotype pattern i.e., star-like topology with a central most-frequent haplotype. However, statistical comparisons of molecular genetic indices: haplotype diversity (H), nucleotide diversity (π), number of polymorphic sites (S), and average number of nucleotide differences (K) showed significant differences between the studied species. The population of E. aenea was significantly more diverse than L. perrisi (P < 0.05, Wilcoxon signed-rank test, Fig. 4).

The W Carpathian population of E. aenea shares haplotypes with locations in Romania, Bulgaria, Finland, Germany and France. 13 COI haplotypes of E. aenea were identified within 276 individuals collected from 47 localities in the W Carpathians (Fig. 2A). The haplotype diversity was 0.336. Considerable genetic homogeneity of the E. aenea population in the W Carpathians resulted from the wide distribution of the dominant haplotype Ea1. The haplotype map (Fig. 2B) revealed that the majority of haplotypes present in southern part of the Carpathians Arc (Romania) and the haplotypes of the Balkan region (Bulgaria) were not recorded from the W Carpathians. The exception was Ea14 shared between one stream in Bulgaria and one spring (V050) in Slovakia (SOM). Individuals from Germany shared the haplotype Ea3 with a single locality in the geomorphological unit FTA (V070). In addition to dominant haplotype Ea1, another five haplotypes were found in FTA and seven in SOM. Haplotypes Ea5 and Ea11 were private for W Carpathians and each occurred in one spring of FTA (V009, V086). The private haplotypes of SOM included Ea7, Ea12 and Ea13, while all of them were located in the springs of SOM1 subunit (V038, V048). Besides that, one spring of SOM1 (V043) shared haplotype Ea8 with the spring of the geomorphological unit FTA (V020). In geomorphological unit WB, four haplotypes were found, while Ea9 was found only at two localities (CZ03, CZ05). On the contrary, the haplotype Ea4 was common in the SOM, FTA, CB and in Romania (Figs. 2A, 2B).

The W Carpathians population of L. perrisi was genetically more homogeneous. Eight haplotypes with a haplotype diversity of 0.007 were found at 43 localities (245 sequences, Fig. 3A). A group of five haplotypes (Lp5, Lp6, Lp7, Lp8, Lp9) recorded in Romania and Bulgaria was highly divergent from the group found in the W Carpathians (Fig. 3B). Haplotype Lp1 dominated in all geomorphological units of the W Carpathians; all German sequences also belonged to this haplotype. The presence of private haplotypes was lower compared to E. aenea: Lp14 from the one spring of SOM (V088), Lp10 in a stream of WB (CZ06) and Lp11 from one spring of geomorphological unit FTA (V070). Besides that, two more haplotypes (Lp3, Lp13) were present in the FTA. Lp13 was shared with the locality of the different geomorphological unit SMC (SEL1). Lp3 was also present in unit WB, in addition to FTA (BEL3). The haplotype Lp2 was detected in a stream (CZ01) of WB and occurred also in one stream (KRV1) located in VM (Figs. 3A, 3B).

The AMOVA showed that most of the observed molecular variance in the W Carpathian populations of both elmid species is generated predominantly within subpopulations (single localities). However, in E. aenea, the molecular variation among subpopulations within geomorphological subunits is more than twice (34.99%) compared to L. perrisi (12.86%) (Table 1). The same results were provided by the AMOVA according to the river basins, where there is also almost no variation associated with the above sub-population (locality) level (Fig. S1, Table S2).

The F_ST values indicate different levels of genetic differentiation between the E. aenea localities 0–0.8 (Fig. 5A) compared to L. perrisi with 0–0.38 (Fig. 5B). The F_ST values of E. aenea suggest that springs V009 (FTA1 - Little Carpathians), V038, and V048 (SOM1 - Slovak Karst) have relatively high pairwise differences in allele frequency, but some level of genetic connectivity cannot be refused. Pairwise comparisons of differences in the frequency of alleles that include mentioned springs are largely significant (P < 0.05). On the other hand, none of L. perrisi F_ST values are statistically significant.

Tests of isolation by distance between springs of both species revealed a positive correlation (Mantel test: E. aenea - r = 0.313, P = 0.000; L. perrisi - r = 0.4122, P = 0.039). Although, only marginally positive but statistically significant correlations suggest a slight structuring effect of the geographical distance among springs in both species. Additionally in case of E. aenea, the spatial autocorrelation was also significant (P = 0.0009, Fig. S2), for L. perrisi it was impossible to calculate it due to lack of data. The genetic distance of E. aenea and L. perrisi in streams was not significantly correlated with the geographic distance (Mantel test: r =  − 0.071, P = 0.153; r = 0.0574, P = 0.115).

Both species were characterized by the statistically significant, negative Fu’s Fs, Tajima’s D and Fu and Li’s D neutrality test values (Table 2). This indicates a recent change in population size of both species. The mismatch distribution analysis suggested a population expansion event for both species, which was indicated by the unimodal shape of the mismatch distribution plot, a small SSD value, and a non-significant p-value (Fig. 6). The eBSP showed a signal of population growth in both species, although the time and character was different. The W Carpathians population of E. aenea (Fig. 7A) started to expand demographically roughly ca. 3,000–2,500 years ago, whereas the population of L. perrisi expanded relatively sharply around 8,000 years ago (Fig. 7B).

Discussion

This study was focused on the two oligo-stenotherm riffle beetles, Elmis aenea and Limnius perrisi (Elmidae). The E. aenea occurred predominantly in karst springs and was rarely found in streams of the W Carpathians, while the distributional pattern of L. perrisi was opposite. This contradicts the previous claims about their common occurrence and similar biotope preference (Moog & Jäch, 1995; García-Criado, Fernández-Aláez & Fernández-Aláez, 1999). Differences in distribution probably can be also explained by altitude, flow type or different ecological demands (Illies & Botosaneanu, 1963). However, according to several studies, E. aenea is more sensitive to harsher conditions resulting from changes of the aquatic environment, manifested, for example, by the loss of macrophytes and moss (Maitland, 1967; Bradley & Ormerod, 2001; Hoffsten, 2003). These findings may explain much greater affinity of E. aenea to springs that generally, with respect to chemical, physical and trophic conditions, are more stable ecosystems compared to other lotic habitats (Minshall & Winger, 1968; Odum, 1971; Butler & Hobbs, 1982; Cushing & Wolf, 1984; Glazier & Gooch, 1987; Gooch & Glazier, 1991). This suggests that karst springs ensured a suitable environment for survival of some aquatic species even during the ice age (Thorup & Lindegaard, 1977). It supports the dinodal hypothesis (Malicky, 1983; Malicky, 2000) proposing that suitable aquatic habitats, persisted throughout the Pleistocene within the periglacial area (dinodal biome), providing suitable conditions for the survival of specialized oligo-stenotherm communities in Central Europe. However, based on our data from the W Carpathians only, the dinodal hypothesis cannot be unequivocally confirmed or refuted, but it clearly opens up new questions in the field of historical-molecular patterns of elmid species in the W Carpathians.

Different results of the Bayesian Skyline Plot analyses between E. aenea and L. perrisi confirmed an exceptional position of the springs. The springs could have a special status in terms of providing stable environmental conditions irrespective of the climatic changes even during the glacial and interglacial periods which did not provoke a dramatic decline or increase of the E. aenea population size in the W Carpathians. In contrast, the populations of L. perrisi, occurring predominantly in streams, began to expand rapidly after the LGP. At the beginning of the Holocene (about 11.5–7.5 ka), a thermal maximum was recorded, which probably enhanced the expansion of species (Dabkowski et al., 2019), which corresponds to sudden expansion of L. perrisi. At that time, local W Carpathian glaciers disappeared completely (Lindner et al., 2003), which led to opening of new migration routes and likely also accelerated species dispersal. Early-Holocene warming is thought to be a major driving force for population divergence in temperate species (Hewitt, 1999). On the other hand, differences in genetic diversity among species, as recorded between E. aenea and L. perrisi may be also influenced by variation in diversification rates (Ricklefs, 2007; Stadler, 2011).

The E. aenea population has occurred in springs at a significantly higher rate, corresponding to its higher molecular diversity compared to L. perrisi that prefers streams and its population is much more uniform. In line with our results, populations of two cofamilial caddisfly species in south-eastern UK showed contrasting genetic patterns. Polycentropus flavomaculatus showed much more pronounced genetic structure than Plectrocnemia conspersa in the same region (Wilcock et al., 2007). In another study on caddisflies of the Central European highlands, Drusus discolor contained three times more haplotypes than Hydropsyche tenuis. Such findings suggest that the isolation of D. discolor populations in Central Europe is stronger and persists for a longer time than in H. tenuis (Lehrian, Pauls & Haase, 2009). Both cases, similarly with species studied herein, confirm that related and co-occurring species may currently have significantly different patterns of molecular diversity, reflecting the different phylogeographical histories of the species and their different autecological traits (Wilcock et al., 2007; Lehrian, Pauls & Haase, 2009).

Compared to aquatic species occurring in streams, the species preferring springs are generally unable to spread extensively and likely persisted at the foothills of mountains during unfavorable climatic conditions (Schmitt, 2007). As a consequence, many of the geomorphological units of the European mountain systems have their own genetic lineages or at least private haplotypes. In our study, we recorded significantly higher values of molecular diversity and higher number of private haplotypes in E. aenea. In addition, analysis of spatial autocorrelation for E. aenea in springs was significant and consistent with these results. This suggests subpopulations in springs persisted in the study area for a longer time, and are relatively isolated. Conversely, stream subpopulations are more homogeneous and smaller, suggesting that they are probably more recent and are being re-created when environmental conditions improve. Valuable examples documenting the importance of the W Carpathians in terms of biodiversity richness were recent discoveries of local endemism of cold-adapted gammarids from Gammarus balcanicus (Mamos et al., 2014; Mamos et al., 2016) and Gammarus fossarum species complexes (Copilaş-Ciocianu et al., 2017) or caddisfly species Drusus discolor. The latter persisted in the Tatra Mts in numerous refugia over multiple glacial cycles, allowing many local endemic clades to form (Pauls, Lumbsch & Haase, 2006). In the case of E. aenea, the two localities in the Slovak Ore Mts (SOM: V038, V048) and one locality in the Fatra-Tatra area (FTA: V009) had remarkably high F_ST values, suggesting that some W Carpathian springs could constitute Pleistocene refugia. According to several studies, the role of the W Carpathians as a glacial refugium (Jamřichová, Potučková & Horsák, 2014; Mráz & Ronikier, 2016; Jamřichová, Petr & Jiménez-Alfaro, 2017) for various species or genetic lineages is undoubted (Pinceel et al., 2005; Magri et al., 2006; Wielstra, Babik & Arntzen, 2015; Mamos et al., 2016; Copilaş-Ciocianu et al., 2017). However, to test whether this hypothesis also applies to Elmidae riffle beetles requires further study in a broader geographical context.

Overall, the genetic differences between populations from different geomorphological subunits of the W Carpathians were very low in both elmid species. However, higher F_ST values in E. aenea correlated with the results of AMOVA. Genetic differentiation among E. aenea subpopulations within geomorphological subunits was relatively high (34.99%). This indicates that there are some well-pronounced differences in genetic composition among most of the spring subpopulations of E. aenea within each geographical unit. Similar results emerged from the study on the black fly Prosimulium neomacropyga in the US Southern Rockies ecoregion with alpine tundra streams, where the differences among streams within the region were 24.58% (Finn et al., 2006). In both elmid species, only a single haplotype was abundant and widespread along the W Carpathians, surrounded by several rare peripheral haplotypes in a star-shaped topology. Similarly, lack of deeper genetic population structure was also found in the W Carpathian populations of the blackfly Simulium degrangei (Jedlička et al., 2012). The maintenance of intraspecific genetic diversity is generally very important for the adaptation potential and long-term survival of species (Spielman, Brook & Frankham, 2004; Frankham, 2005). However, prolonged persistence is possible even despite low levels of genetic diversity (Johnson et al., 2009). Relatively homogeneous population patterns of both studied riffle beetles may reflect their short history in the W Carpathians. The comparatively low genetic differentiation among populations of trickle midges (Diptera: Thaumaleidae) in Northern Europe was also explained by relatively recent, possibly post-glacial dispersal (Haubrock et al., 2017).

Conclusions

In conclusion, it seems that different habitat preferences of the two related aquatic beetle species E. aenea and L. perrisi preserved their similar population-geographical patterns, but shifted their molecular diversity, as well as the time and character of their distribution in the W Carpathians. E. aenea, with higher molecular diversity, occurred mainly in the springs compared to the genetically more homogenous population of L. perrisi that was found mostly in streams. These findings support the attribution of the W Carpathian springs to potential refugia with a suitable environment allowing for survival of aquatic biota even during the unfavorable climatic conditions through geological ages and maintaining or even developing its intraspecific genetic diversity. In addition, an isolation of the W Carpathians springs is also indicated by the significant results of Mantel and spatial autocorrelation analysis in E. aenea. This study added new information about understudied riffle beetle fauna of one of the world’s biodiversity hotspots, the W Carpathians. However, further studies should include more samples from Southern and Eastern Europe in order to understand the holistic biogeographic pattern of the target species and spring fauna in general.

Supplemental Information