The karst areas in the Dinaric region of the Western Balkan Peninsula are a hotspot of freshwater biodiversity. Many investigators have examined diversification of the subterranean freshwater fauna in these karst systems. However, diversification of surface-water fauna remains largely unexplored. We assessed local and regional diversification of surface-water species in karst systems and asked whether patterns of population differentiation could be explained by dispersal—diversification processes or allopatric diversification following karstrelated microscale vicariance. We analyzed mitochondrial cytochrome c oxidase subunit I (mtCOI) sequence data of 4 caddisfly species (genus Drusus) in a phylogeographic framework to assess local and regional population genetic structure and Pliocene/Pleistocene history. We used BEAST software to assess the timing of intraspecific diversification of the target species. We compared climate envelopes of the study species and projected climatically suitable areas during the last glacial maximum (LGM) to assess differences in the species climatic niches and infer potential LGM refugia. The haplotype distribution of the 4 species (324 individuals from 32 populations) was characterized by strong genetic differentiation with few haplotypes shared among populations (16%) and deep divergence among populations of the 3 endemic species, even at local scales. Divergence among local populations of endemics often exceeded divergence among regional and continental clades of the widespread D. discolor. Major divergences among regional populations dated to 2.0 to 0.5 Mya. Species distribution model projections and genetic structure suggest that the endemic species persisted in situ and diversified locally throughout multiple Pleistocene climate cycles. The pattern for D. discolor was different and consistent with multiple invasions into the region. Patterns of population genetic structure and diversification were similar for the 3 regional endemic Drusus species and consistent with microscale vicariance after the onset of intensified karstification in the Dinaric region. Karstification may induce microscale vicariance of running surface-water habitats and probably promotes allopatric fragmentation of stream insects at small spatialscales.
Vicariance and dispersal abilities are important historical determinants of a species’ geographical range (Beebee and Rowe 2004). Isolation of populations caused by the emergence of geographical or environmental barriers has resulted in disjunct distribution patterns in many taxa (e.g., Schneeweiss and Schönswetter 2010, Viruel et al. 2011). Thus, vicariance generally is considered an important process in many allopatric speciation events, mostly on larger geographic scales (e.g., Tethyan vicariance: Hrbek and Meyer 2003, Bunje and Lindberg 2007). Dispersal across pre-existing barriers, e.g., range expansions into new, isolated habitat, also can mediate speciation (Platnick and Nelson 1978). At larger scales, past climate change has influenced dispersal and range shifts (Schmitt 2007) and can promote dispersal-mediated diversification in various taxa (e.g., Erkens et al. 2007, Milá et al. 2007). Differentiation of species and populations on small scales also are often associated with limited dispersal (e.g., Brändie et al. 2007, Johannesen et al. 2010), and small-scale migrations have frequently resulted in fragmentation or diversification of populations (e.g., Schmitt et al. 2006, Margraf et al. 2007).
The general importance of small-scale vicariant and dispersal-mediated speciation in areas of high regional endemism is understudied. We define small-scale vicariance as changes in the landscape that lead to fragmentation of suitable habitats and, thus, populations at the local and regional scale. Dispersal-mediated speciation is also an allopatric diversification process but is initiated by regional dispersal into newly colonized, isolated habitats. Both processes play varying roles in different aquatic and terrestrial taxa. In a lineage of aquatic beetles with many narrowly distributed allopatric endemics in southwestern Europe, the 2 processes influenced diversification by varying degrees in different clades (Ribera et al. 2011). Both processes also were identified as having shaped distribution patterns of stygobionts, i.e., obligate subterranean-dwelling aquatic species, in diverse karst regions of the world. The relative importance of the 2 processes in karst areas differs from region to region (Culver et al. 2009). Differences in the geological history of particular regions generate varying opportunities for vicariance, and differences in life-history traits of terrestrial and aquatic subterranean fauna also are important factors governing the role of the 2 processes (Culver et al. 2009).
Karst areas generally have high degrees of endemism, particularly of freshwater and subterranean fauna (e.g., Gibert and Culver 2009, Freyhof and Brooks 2011). The karst area in the Dinaric region of the Western Balkans is the world's richest region of subterranean biodiversity, especially the stygobiont fauna (Gottstein-Matočec et al. 2002, Gibert and Culver 2009). Karstification (dissolution of soluble rocks, e.g., limestone and dolomite) generally changes topography and hydrology in karst areas (Bosák 2008) and leads to a reduction of surface waters and high levels of local topographic relief with deeply incised canyons. For the subterranean fauna, karst development and increases in cave density create increased opportunities for migration (Culver et al. 2009). In contrast, for organisms inhabiting surface waters, karstification increases habitat fragmentation. Flight appears to be the major dispersal mechanism in stream insects, except in cases where adult flight capability is very poor or where terrain is extremely mountainous (Hughes et al. 2008). Thus, reduction of surface waters and parallel development of rugged and steep topography could limit both habitat availability and adult dispersal among streams. Karstification may drive small-scale vicariance of stream insect populations, thereby leading to local allopatry (e.g., Klobučar et al. 2013). Thus, local vicariance rather than dispersal may drive diversification and explain the high number of microendemic surface-water species in karst regions.
In particular, the Dinaric region of the Western Balkans has high diversity of surface-water invertebrate species with small ranges (Schmidt-Kloiber and Hering 2012). Intensive morphological and molecular genetic analyses (Marinlcović-Gospodnetić 1971, 1976, Malicky 2005, Previsie et al. 2009, Kucinic et al. 2011) show that many endemic caddisfly species with very localized distribution ranges exist in the region. Karst development at the landscape scale has been studied intensely and is quite well known in the Dinaric region (Mihevc et al. 2010) because this region of the Western Balkans is the locus typicus of the global karst. Thus, microendemic caddisfly diversification in this region is a good model for testing hypotheses on small-scale vicariance in freshwater organisms, how it may be linked to karst development, and its influence on surface-water availability.
Comparative phylogeographic studies of codistributed species have provided insights into the evolutionary history of lineages and species in many regions, e.g., by inferring common biogeographic breaks, recolonization routes, or locations of refugia (e.g., Avise 1992, Schmitt 2007, Lehrian et al. 2009). However, phylogeography also can be used to assess the historical responses to environmental change of closely related taxa. In particular, such studies can be used to examine evolutionary constraints upon which intraspecific processes act in shaping the genetic structure of populations (Hodges et al. 2007, Pauls et al. 2009).
Comparative phylogeographic studies are lacking in the Western Balkan region, and common phylogeographic patterns of biota inhabiting this region have not yet been documented. We used 4 caddisfly species of the subfamily Drusinae (Trichoptera:Limnephilidae) that are closely related and partially codistributed (Fig. 1A) to examine spatial and temporal effects of small-scale vicariance at both intra- and interspecific levels. To date, 95 Drusinae species have been described. Some of them are widely distributed, but most are endemics limited to a single or few mountain ranges across the Eurasian mountains (Malicky 2005, Graf et al. 2008, Oláh 2010, 2011, Kučinić et al. 2011). The Pyrenees, the Alps, the Dinaric region of the Western Balkans, the Carpathians, and eastern Mediterranean are particularly diverse regions (Malicky 2005). Of the 21 Drusinae species recorded from the Dinaric region of the Western Balkans, 12 are endemic (Marinković-Gospodnetić 1971, 1976, Graf et al. 2008, Kučinić et al. 2011). Thus, the group presents a valuable model for studying evolutionary processes like speciation and diversification. The species studied here are cold-water stenotherms restricted to permanently flowing headwaters and springs, as are most Drusinae (Graf et al. 2008). Dispersal abilities of Drusinae are poorly known. Behavior plays an important role in dispersal (Engelhardt et al. 2011, Müller-Peddinghaus 2011), so assessing their dispersal capacity is difficult. Wing morphology blueprints suggest that Drusinae are relatively poor dispersers (Müller-Peddinghaus 2011). However, Drusinae are caught regularly in light traps, but the distances from the sources and frequency can vary among species (Malicky 1987). Differences in population genetic structure have been attributed to varying levels of dispersal among closely related Drusinae (Pauls et al. 2009).
Figure 1.
Distribution (A), sampling and population genetic structure (B-E) of 4 Drusus species (Drusinae:Limnephilidae: Trichoptera) within the Dinaric region of the Western Balkans. Sampling sites and median joining (MJ) mitochondrial cytochrome c oxidase subunit I (mtCOI) haplotype networks (shown as insets) are shown separately for Drusus discolor (B), Drusus schmidi (C), Drusus croaticus (D), and Drusus krusniki (E). Populations are coded by country (AT: Austria, BIH: Bosnia and Herzegovina, HR: Croatia, KS: Kosovo, MK: Macedonia, MN: Montenegro) and colored by sampling locality. Haplotypes in the MJ network are colored according to geographic origin of specimen(s) with a given haplotype. Haplotype diameter is relative to the frequency of specimens with a given haplotype. Genetic distances between haplotypes are to scale, unless indicated. Distances >1 basepair (bp) are indicated by the number of changes.
We extended previous studies of geographic structure and timing of divergences in Dinaric Drusus species (Previšić et al. 2009, Kučinić et al. 2011) by comparing population genetic structure and the phylogeography of 4 closely related species based on range-wide intraspecific data sets. We assessed whether small-scale vicariance or dispersal were important factors in the diversification of aquatic insects in the Western Balkans. We hypothesized that karstification-mediated fragmentation of suitable surface-water habitats and microscale vicariant diversification were the primary processes driving speciation of highland caddisflies in the Dinaric region of the Western Balkans. We used mitochondrial deoxyribonucleic acid (mtDNA) sequence data to analyze population genetic structure, date intra- and interspecific divergence of our target species and intraspecific lineages, and link these with geological events, particularly the onset of intensified karstification in the late Pliocene and early Pleistocene. We tested an alternative scenario based on climate-change-induced large-scale migrations and species range expansion/regression cycles that could have led to population fragmentation and diversification in the region. We used species distribution models (SDM) to test this scenario by assessing the suitability of the last glacial maximum (LGM) climate conditions for each of the target species in the Western Balkans. We predicted that the target species found suitable climatic conditions in their current ranges during the LGM and were not subject to large-scale migrations and species range expansion/ regression cycles. If the species were subject to large-scale migrations, we would expect patterns indicative of secondary contact or admixture, particularly considering the close proximity of many present-day populations. Our study offers more detailed insights into the validity and generality of the refugia-within-refugia hypothesis in the Western Balkans (Previšić et al. 2009) and contributes to our understanding of how microscale vicariance and dispersal-mediated diversification affect diversity in coldadapted highland aquatic insects.
METHODS
Study species
Drusus discolor (Rambur, 1842) has a disjunct distribution across most of the European mountain ranges, including the Dinaric Alps (Pauls et al. 2006, 2009) (Fig. 1B) and is generally restricted to headwaters at higher altitudes. Genetic differentiation and numerous regionally endemic lineages were observed in D. discolor across Europe (Pauls et al. 2006, 2009). Populations from the Western Balkans were previously not available for study, and we lack an understanding of population genetic structure in this region. Drusus schmidi Botosaneanu, 1960 is endemic to the Dinaric region of the Western Balkans and also inhabits headwater streams (Graf et al. 2008) (Fig. 1C). The microendemic Drusus croaticus Marinkovic-Gospodnetic, 1971 (Fig. 1D) and Drusus krusniki Malicky, 1981 (Fig. 1E) have very restricted distribution ranges in the northwestern and eastern part of the Dinaric region of the Western Balkans, respectively (Marinković-Gospodnetić 1971, Kući-nić et al. 2008, Oláh 2010, 2011). These 2 species inhabit cold springs and show no specialization to certain spring types (AP, WG, MK, personal observation). A previous study of the population genetic structure of D. croaticus revealed very high genetic differentiation on a small geographical scale (Previšić et al. 2009).
Sample collection, DNA extractions, polymerase chain reaction amplifications
We collected specimens of the target species from stream substrates and riparian vegetation using hand nets. We extracted DNA from specimens using the DNeasy Tissue Kit or the QIAmp Micro Kit (both Qiagen, Hilden, Germany). We amplified a 541-basepair (bp) fragment of mtDNA from the cytochrome c oxidase subunit I (mtCOI) gene region. Details of molecular methods are given in Pauls et al. (2006) and Previšić et al. (2009).
Population genetic structure
We estimated population genetic structure using analysis of molecular variance (AMOVA; Arlequin, version 3.11; Excoffier et al. 1992), which partitions genetic variance among and within populations. AMOVA was calculated based on Kimura-2-parameter (K2P) distance and 10,000 bootstrap replicates. We also calculated exact tests of population differentiation (Raymond and Rousset 1995) with 100,000 steps in the Markov chain and a burn-in of 10,000 steps, as well as pairwise FST based on K2P distance and 10,000 bootstrap replicates. Both analyses were run in Arlequin (version 3.11; Excoffier et al. 2005). We calculated Median Joining (MJ) networks (Bandelt et al. 1999) using default settings in Network (version 4.5.0.0; Fluxus Technology Ltd., Suffolk, UK) to illustrate diversity and geographic distribution of haplotypes for all species. Median Joining Networks were drawn in Network Publisher (version 1.3; Fluxus Technology Ltd.), and the figures prepared in Adobe Illustrator (version 15.0; Adobe, San Jose, California).
Estimation of divergence time
For estimates of divergence time we included 15 endemic Western Balkan Drusus species of monophyletic origin (AP, SUP, unpublished data). We estimated interspecific divergence times among these 15 species and among intraspecific lineages of D. croaticus, D. discolor, D. krusniki, and D. schmidi from mtCOI sequences with a Bayesian Monte Carlo Markov Chain (MCMC) coalescent method implemented in BEAST (version 1.6.1; Drummond and Rambaut 2007). For all data sets, Hasegawa-Kishino-Yano + Invariant (HKY+I) was identified as the best-fitting model of DNA substitution using the Akaike Information Criterion (AIC) test implemented in MrModeltest (version 2.2; Nylander 2004). Reliable calibration points are not available for these caddisflies, so we used the sequence divergence rate for Lepidoptera, the closest relative to modern Trichoptera, of 2.3% per million years for mtCOI sequences (Brower 1994). We set the tree prior to a Yule process and a constant-size coalescent prior for analysis of all 15 Western Balkan Drusus species and populations of the 4 target species, respectively. We ran the analyses with an uncorrelated relaxed molecular clock with 2 independent chains of 50 million generations, sampling every 1000 generations, and a burn-in of 5 million generations. Convergence of independent chains and adequate mixing of all parameters were checked in Tracer using the Effective Sample Size diagnostic (version 1.4; Drummond and Rambaut 2007). After the independent runs were combined in the program LogCombiner (version 1.6.1; Drummond and Rambaut 2007), we used TreeAnnotator (version 1.6.1; Drummond and Rambaut 2007) to build the maximum clade credibility tree.
Species distribution modeling
We took distribution data for all species from our own collections and from published data based on expert identifications ( Table S1 (TableS1.pdf)). The number of unique localities was 13 for D. schmidi, 17 for D. krusniki, and 19 for D. croaticus. Drusus discolor was modeled based on its entire European range (528 sites), including 17 sites from the target region ( Table S1 (TableS1.pdf)). We did species distribution modeling (SDM) for each of the 4 target species in MaxEnt (version 3.3.3; Phillips et al. 2006) using environmental data obtained from the WorldClim database ( http://www.worldclim.org/) at 5-km2 spatial resolution. MaxEnt outperforms most other modeling applications when using presence-only data (Elith et al. 2006) and is least affected by small sample sizes (Pearson et al. 2007, Wisz et al. 2008) and georeferencing errors (Graham et al. 2008). We produced several models based on different sets of climatic parameters: 1) all 19 bioclimatic variables, 2) selected variables based on expert knowledge on the ecology of the focal species, excluding highly correlated variables (r > 0.7), and 3) selected variables based on model importance of the initial (full) model (again avoiding high correlations). A full list of variables for each modeling run is provided in Table S2 (TableS2.pdf). For each species, we ran 10 replicates with cross-validation. To assess potential range shifts since the LGM, we projected the climatic envelopes onto 2 LGM scenarios (Community Climate System Model [CCSM] and Model for Interdisciplinary Research of Climate [MIROC], both available from the WorldClim database). Models for wide-spread D. discolor were built over its entire range (Western Europe), but only the Western Balkan region is shown here for comparison with the other taxa.
Niche overlap between species (at the present and for each LGM scenario) and temporal changes for each species (i.e., overlap between present-day and LGM patterns) were assessed using 3 different similarity statistics: Schoener's D (Schoener 1968), the I statistic (Warren et al. 2008), and the relative rank test (Warren and Seifert 2011). Schoener's D and I calculate the difference in standardized suitability scores for each grid cell, whereas the relative rank test provides an estimate of the congruence of relative ranks of suitability for each grid cell. All 3 measures range between 0 (nonoverlapping niches) and 1 (identical niches).
RESULTS
We sequenced 541 bp of mtCOI from 324 specimens of D. croaticus (Nind = 116, Npop = 11), D. discolor (Nind = 87, Npop = 11), D. krusniki (Nind = 88, Npop = 6), and D. schmidi (Nind = 31, Npop = 5) (Fig. 1B-E). We also generated mtCOI sequences from 14 additional Drusinae species endemic to the Western Balkan region ( Table S1 (TableS1.pdf); Gen-Bank Accession Numbers FJ002601-FJ002721, GQ470606, GQ470607, and KC881310-KC881525). Some sequences used in our study were published by Previšić et al. 2009 and Kučinić et al. 2011 (see Table S1 (TableS1.pdf) for details).
All populations showed significant population genetic structure (ΠST = 0.912–0.959, all p<0.0001; Table 1). Strong population genetic structure also was confirmed by population pairwise estimates of FST and exact tests of population differentiation ( Table S3 (TableS3.pdf)).
All species showed deeply divergent, geographically restricted clades (Fig. 1B-E). All haplotypes in all 4 species were regionally endemic, though some (16%) were shared among neighboring populations. Clades were minimally separated by 2 bp changes in D. krusniki, 6 bp changes in D. schmidi, and 8 bp changes in D. discolor and D. croaticus (Fig. 1B-E). The most divergent haplotypes within our target species were separated by 20 bp changes in D. krusniki, followed by 14 bp changes in D. croaticus. In D. discolor, 3 regionally differentiated clades were associated with the Karawanken Mountains; the northwestern Dinaric Alps; and the southeastern Dinaric Alps and Pelister and Kožuf Mountains in Macedonia (Fig. 1B). The clades showed shared haplotypes populations, but no shared haplotypes among these areas. The 3 clades did not form a monophyletic lineage among all known mtCOI haplotypes of D. discolor from across Europe (SUP, unpublished data). The other species were more narrowly distributed (Fig. 1C-E), but they exhibited much deeper divergence among geographically proximate populations, particularly D. krusniki and D. croaticus (Fig. 1D, E).
Table 1.
Results of the analysis of molecular variance (AMOVA) of 4 Drusus species within the Dinaric region of the Western Balkans.
Figure 2 shows the majority-rule consensus tree of 2 independent BEAST runs for the individual target species and a group of the 15 Western Balkan endemics. Analysis of the log files from independent BEAST runs indicated parameter convergence and good chain mixing as measured by effective sample size (ESS) values for all model parameters (data available online from the BiK-F MetaCat Database from http://dataportal-senclcenberg.de/bikfdata/). Based on the ultrametric trees, the population divergence of all 4 study species was estimated to fall within the Pleistocene, prior to the LGM, but the confidence intervals for D. krusniki and D. croaticus date back to the late Pliocene (Fig. 2). In general, the populations of the narrow-range endemics D. croaticus and D. krusniki showed older divergences than D. schmidi or D. discolor (Fig. 2). The monophyletic clade of the 15 Western Balkan endemics was estimated to have diverged into 2 groups during the late Miocene. Divergences between closely related species within each of these 2 groups were estimated to range from early Pliocene to late Pleistocene, with most splits dating to the Pleistocene (Fig. 2).
Distribution modeling under current climatic conditions provided patterns of habitat suitability largely concordant with known species distributions throughout the region ( Fig. S1 (FigureS1.pdf)). Model performance, evaluated using the area under the curve (AUC) value derived from the receiver operating characteristic (ROC) curve was high for all species and across models. All mean AUC values for the replicate runs were >0.938. Only a few climatically suitable habitats were identified outside the known species ranges. Projected LGM distributions of suitable climatic conditions largely overlapped with present-day patterns for most species ( Table S4 (TableS4.pdf)) and generally indicated suitable ranges during the LGM in or close to the species’ current ranges. In most cases, small differences exist between the 2 LGM scenarios, but the differences are not consistent in the sense that one or the other scenario always infers wider/narrower ranges. In D. croaticus, the CCSM projections showed a strong deviation from the current and MIROC LGM projections of suitable conditions. No suitable conditions were projected in the northwestern part of the study region based on CCSM.
The climatic envelopes of the 4 species varied distinctly (Fig. 3). Drusus discolor occupied the greatest breadth of temperature and precipitation conditions. The occurrences of D. schmidi and D. krusniki were associated with a relatively broad mean annual temperature range (∼3– 10°C), but a narrow annual precipitation range (∼1000– 1200 mm). Drusus croaticus exhibited the narrowest climatic envelope characterized by high annual precipitation (∼ 1200–1400 mm) and relatively warm mean annual temperature (∼7–ll°C).
Niche overlap between different Drusus species within the Western Balkans indicates taxon-specific changes since the LGM. The overlap between D. krusniki and D. schmidi remained relatively constant, but both increases and decreases can be seen for the other species ( Table S5 (TableS5.pdf)). Since the LGM, the overlap between D. croaticus and all other species increased, whereas the overlap of both D. krusniki and D. schmidi with D. discolor decreased. The degree of pairwise range overlap differed between the 2 LGM scenarios and the niche overlap statistics used. Nevertheless, the overall patterns of change in niche overlap through time described above remained stable regardless of the scenario or the statistic.
DISCUSSION
Is diversification driven by karstification or climate-change-induced dispersal-migration?
Both karstification and Pleistocene climate change have been inferred as important processes of diversification in the Dinaric region of the Western Balkans (e.g., Sket 1999, Previšić et al. 2009, DeFaveri et al. 2012). The patterns of diversification we observed in our study in widespread and endemic species and at the inter- and intraspecific level are generally in concordance with the onset of the Pleistocene, i.e., major climate change and intensification of the karstification process in the region (Sket 1999, Mihevc et al. 2010). Some of our estimates of divergence times differ from previous estimates of divergence of Drusus endemics in the Dinaric region of the Western Balkans. Previous dates were estimated from a more limited data set and with different models applied for divergence estimation (Previšić et al. 2009). Both factors can affect divergence time estimates (Drummond et al. 2006). How-ever, results of previous studies and our results support regional speciation and differentiation beginning in the late Pliocene through the Pleistocene. In general, divergence estimates must be interpreted with care. Different taxa can exhibit different evolutionary rates at the same locus, and differences often are higher within than among species (e.g., Papadopoulou et al. 2010). The rate we applied here is from a quite divergent insect group and could be slower than the true rate. Therefore, we interpret the data cautiously and with attention to the large error ranges indicated in the BEAST estimates of divergence time. Basal intraspecific divergence times for all 4 species are within the late Pliocene to mid-Pleistocene range and would remain so, even if we used a much higher rate like the one proposed by Papadopoulou et al. (2010). On the other hand, the shallowest intraspecific splits could shift to the Pleistocene/Holocene boundary.
Figure 2.
Estimates of divergence times for 15 endemic Western Balkan Drusus species of monophyletic origin and intraspecific divergence of the 4 target species based on BEAST analysis. Individual haplotypes were collapsed by clades. Clades are coded according to their geographic origin (see Fig. 1). Green bars indicate 95% highest posterior density (HPD) intervals for divergence times of individual nodes (Node *: mean age estimate = 5.288 Mya and 95% highest posterior density [HPD] interval = 7.77– 3.31 Mya; Node **: 95% HPD interval = 3.72–1.61 Mya). Vertical grey lines show 250 ky increments. Epoch boundaries are shown according to Litt et al. 2008 (L.P. = late Pleistocene, and H. = Holocene).
Figure 3.
Scatterplot of temperature and precipitation associated with the occurrences of the 4 study species. Points represent collection localities of each species, histograms show the frequency distribution of each species with regard to annual precipitation (y-axis) and mean annual temperature (х-axis) derived from the collection localities.
In general, distribution patterns and genetic differentiation in many taxa have been attributed to the interplay of vicariance and dispersal (e.g., Veith et al. 2003, Ribera et al. 2011), but information on the importance of these processes in taxa inhabiting the Balkans is still limited (but see Culver et al. 2009, Triponez et al. 2011). All our study species exhibit a highly differentiated population genetic structure that supports a scenario in which allopatrically isolated lineages formed by vicariance on a very small geographic scale. Highly differentiated lineages on small spatial scales already have been observed in various taxa in the Western Balkan region (e.g., various freshwater fishes: Šlechtová et al. 2004, DeFaveri et al. 2012; crayfish: Klobučar et al. 2013; cave-dwelling bivalve: Bilandžija et al. 2013).
Nevertheless, the isolating effects of different processes are not easy to distinguish. Both vicariance through Pliocene/Pleistocene karstification on the one hand and Pliocene/Pleistocene range shifts in combination with population expansion/regression cycles on the other hand could have resulted in population fragmentation and diversification (e.g., Sanmartín 2003, Triponez et al. 2011). Vicariance associated with karstification could be the case under climatic niche conservatism, which can promote allopatric isolation of lineages that track their climaticniches and can, in some cases, lead to speciation (Kozak and Wiens 2006). On the other hand, under a range regression/expansion scenario with large-scale migration and repeated dispersal events, we would expect to see patterns of secondary contact, particularly at the regional and small scales as have been observed in other regions for D. discolor (Pauls et al. 2006, 2009), other montane caddisflies (Engelhardt et al. 2011), montane mayflies (Theissinger et al. 2011), and crustaceans in the Western Balkans (Zakšek et al. 2009, Wielstra et al. 2010). This pattern is not consistent with any of our study species, where only genetically close haplotypes (<4 bp divergence) are shared among geographically proximate populations. Further-more, the endemic species (D. croaticus in particular) show locally isolated lineages that appear to have been separated for long periods of time (600–300 kya) and almost certainly are older than the last interglacial/glacial cycle (Eem/ Würm cycle in the Alps, 126 to 11 kya; Litt et al. 2008).
In line with the genetic data suggesting long-term persistence of individual populations in isolation, the SDM projections generally support the persistence of all 4 target species within the Dinaric region during the Pleistocene cold stages, without admixture events that disrupted the prevailing population genetic structure. Discrepancies between the 2 LGM scenarios (CCSM and MIROC) are observed in all target species, but are particularly pronounced in D. croaticus. A similar phenomenon has been observed in other studies of aquatic (e.g., Rödder and Dambach 2010) and terrestrial (e.g., Habel et al. 2010, Vega et al. 2010) organisms covering the Dinaric region of the Western Balkans. The differences between the 2 LGM scenarios can be attributed to strong discrepancies in the temperature conditions under the 2 models in the north-western Dinaric region ( Fig. S2 (FigureS2.pdf)). The temperature conditions of the MIROC model are much closer to present-day temperatures than those of the CCSM, particularly in the northwestern part of the Dinaric region of the Western Balkans and, thus, precisely in the small present-day distribution range of D. croaticus. Considering the warmer temperature range of D. croaticus compared with the other species, the CCSM model probably projects LGM temperatures that are too cool in its present-day distribution range. The 2 models differ dramatically in temperature projections in this region, highlighting the regional inconsistencies among global models and their limitations when used at regional scales. Unfortunately, regional climate and LGM models do not exist for the study region.
Our results are based only on the matriline of the study species and should be interpreted with care. The results do suggest that the karstification process played an important role in promoting the microscale allopatric diversification of aquatic organisms caused by reduction of the surface drainage, i.e., loss and fragmentation of suitable habitats. This outcome is in contrast to the situation in the cave fauna, which profited from the increase in subterranean water discharge (Sket 1999). The effects of karst development on hydrology in the region were generally very drastic, but varied locally (e.g., Sket 2002, Prelovsek 2010). In addition to the major influence of extensive karstification, decreases in the Adriatic Sea level during glaciations led to lower groundwater levels (Surić et al. 2005), which may have further reduced surface flows in the region and promoted isolation among surface-water habitats for stream organisms. Thus, the interplay of various geological events and processes resulted in microscale vicariance that may have driven diversification in our study species and, more generally, in the freshwater fauna establishing a center of speciation and a biodiversity hotspot (Freyhof and Brooks 2011, DeFaveri et al. 2012).
Phylogeography of 4 co-occurring caddisflies—new and unique patterns
All 4 species have highly differentiated population genetic structure, but differentiation seems to be strongest among the endemic species. These species exhibit minimal sharing of haplotypes among neighboring populations and very high clade divergence. They are obviously poor dispersers, even at local scales. In contrast, D. discolor appears to have colonized the 3 regions (the eastern part of the Southern Alps including Karawanken; northwestern Dinaric Alps; southeastern Dinaric Alps, Pelister, and Kožuf Mountains in Macedonia) independently because the 3 lineages do not form a monophyletic clade (SUP, unpublished data). Most haplotypes are shared among regional populations in D. discolor, a result indicative of greater dispersal activity and gene flow in this species than the regional endemics. Because of disagreement with the SDM projections, we cannot rule out either explanation for D. discolor with data at hand, i.e., persistence within the Dinaric region during the Pleistocene cold stages and range expansion/regression cycles, or even a more complex scenario involving both survival in localized refugia and large-scale migrations.
Dispersal and vicariance are key processes driving allopatric diversification, but ecological diversification also can be important in the speciation process (Elias et al. 2012). Differences in species ecology and life-history traits (e.g., feeding strategies, habitat preferences) also can shape population histories (e.g., Füreder and Pöckl 2007). For instance, in freshwater shrimps, variations in egg size have been linked with differing dispersal abilities, thus affecting population structure, intraspecific divergence, and geographic distribution (Page and Hughes 2007). Ecological trait information regarding feeding, thermal preferences, climate envelopes, and habitat preferences are available for the target species. Drusus discolor is a filtering carnivore, whereas the remaining species are grazers (Graf et al. 2008). Predatory caddisfly species may be more opportunistic than more specialized grazers in their feeding ecology, thus putting them at a competitive advantage in the establishment process (e.g., Malmquist et al. 1991). This aspect may be particularly important in Drusinae because the clade of potentially competing grazers is much more diverse than the other feeding-type clades (Pauls et al. 2008). All 4 target species are coldwater stenotherms, but their habitat preferences differ markedly (Graf et al. 2008). In addition, they clearly differ in terms of their climatic occurrences (Fig. 3). Species occurrences cannot be linked directly to climatic preferences or climatic niches, but they can provide an indication of climate sensitivity. The endemic species exhibit much narrower temperature and precipitation ranges than D. discolor. This narrowness may be an indication of local adaptation and specialization in the endemics, but is much more likely to be an artefact of inferring the climatic occurrence across the area associated with the smaller range size (Wisz et al. 2008). However, the very narrow precipitation ranges of the endemics and the clear temperature differentiation between D. croaticus and the other 2 endemics suggest some degree of differential climate sensitivity. Species living in karst regions are particularly dependent on precipitation patterns to maintain their above-ground freshwater habitats (e.g., Bonacci 2001). In terms of climate niches, the LGM projections suggest that climate conditions generally remained suitable and, thus, presumably more-or-less stable through time. Changes in niche overlap since the LGM reveal different patterns for the Dinaric endemics. Niche overlap has remained more-or-less constant between D. krusniki and D. schmidi, but overlap between D. croaticus and the other 2 species apparently has increased since the LGM. This increase might suggest that the potential for genetic exchange was probably not higher in the past than under current climatic conditions. Whether the differences in phylogeographic history and present population structure result from differences in dispersal behavior, climatic sensitivity, or other ecological traits cannot be answered fully with the data at hand.
Our study reveals 3 distinct patterns for 3 species. This outcome is in line with previous assessments that phylogeographic patterns of montane aquatic insects seem to be species-specific and that general patterns are rare (e.g., Lehrian et al. 2009, Pauls et al. 2009). The pattern of geographic distribution coupled with the genetic differentiation as observed in Drusinae has not been recorded within the Dinaric region of the Western Balkans. Some common phylogeographic patterns, such as the local differentiation centers in different animals (e.g., Kryštufek et al. 2007, Ursenbacher et al. 2008) and plants (e.g., Surina et al. 2011), or common historical breaks and connections in freshwater fishes (summarized by Perea et al. 2010), have been documented in the region. However, the patterns observed in Drusinae are quite unique, showing that vicariance and dispersal operating in such a geologically complex environment may create distinctive evolutionary patterns, resulting in high biodiversity (e.g., Hewitt 2011).
Conclusions
Our study demonstrates that population genetic structures of small-range endemic aquatic insect species show similar patterns probably linked to microscale vicariance in karst landscapes. SDM projections are informative and are a good supplement to phylogeographic studies, but our results show that the SDM models and LGM projections based on global models should be used with caution for specific areas and at regional scales. In our case, the integrative approach provides further support for the refugia within refugia hypothesis in the Balkans. Moreover, our results suggest that vicariance through land-scape alteration, such as karstification, can play an important role in allopatric diversification, even on very small spatial scales.
ACKNOWLEDGEMENTS
We are very grateful to the following colleagues for helping us in the field or for providing specimens: Astrid Schmidt-Kloiber (Vienna); Matej Faller, Iva Mihoci, Aleksandar Popijač, and Krešimir Žganec (Zagreb); Svjetlana Stanić-Koštroman (Mostar); Vladimir Pešić (Podgorica); Vladimir Krpač (Skopje); and János Oláh (Debrecen). We thank Miklos Bálint (Frankfurt) for helpful discussions and providing the altitude and slope layers for the SDM, and Charles T. Wardwell (St Paul) and Jürgen Otte (Frankfurt) for assistance in the laboratory. We also thank Tamás Mikes (Frankfurt) for comments that improved this manuscript. AP gratefully acknowledges a BiK-Net incoming fellowship. This study was financially supported by the research funding programme “LOEWE - Landes-Offensive zur Entwicklung wissenschaftlich-ökonomischer Exzellenz” of Hesse's Ministry of Higher Education, Research, and the Arts; Croatian Ministry of Science, Education and Sports (Project Nos. 119–1193080–1206 and 119–1193080–3076), the Austrian Science Fund (FWF) (Project No. 23687-B17) and the BioFresh EU project-Biodiversity of Freshwater Ecosystems: Status, Trends, Pressures and Conservation Priorities (contract No. 226874).
AP and SUP designed the research; AP, MK, MK, WG, and HI contributed new reagents/analytic tools; AP, JS, and SUP performed the research, analyzed the data, and wrote the paper. All authors edited the drafts and approved the final version of the manuscript.
