We used 2 caddisflies, Drusus discolor and Drusus romanicus, to test explicitly whether closely related species that occupy similar niches and occur in partial sympatry maintain comparable population structure and share a similar population history. We used mitochondrial sequence data to analyze and compare the population structure and the phylogeography of 105 specimens of D. discolor and 74 individuals of D. romanicus collected in southeastern Europe. We examined the relationship between both species with phylogenetic inference and coalescent modeling and used the results to assign larvae to species. We were able to unambiguously assign larvae to species level based on clearly defined association criteria within a phylogenetic analysis of all specimens. The species were closely related and not reciprocally monophyletic in our haplotype phylogeny. One D. romanicus haplotype from the Bucegi Mountains was nested within D. discolor, a result that suggests isolation with migration, introgression, or incomplete lineage sorting between the 2 species. For each species, we examined population genetic structure with median joining networks, analysis of molecular variance (AMOVA), exact tests of population differentiation, and Mantel tests of isolation by distance. We used tests for selective neutrality (Tajima's D, Fu's F) to infer potential population growth and expansion. Species differed in their genetic population structure. Drusus discolor had haplotype overlap among several mountain ranges in the study region. No D. romanicus haplotypes were shared among any regions examined, and levels of divergence between haplotype clades exceeded those of D. discolor by a factor of up to 2.1. The different degree of population differentiation and divergence of both species probably reflects different Pliocene/Pleistocene population histories and might be related to differences in dispersal capabilities or competitive exclusion of D. romanicus by D. discolor in the mountain ranges north and west of the Western Carpathians. Based on our results, we discuss the importance of the Carpathian Mountains and Bulgarian highlands as Pliocene/Pleistocene refugia and centers of diversification.
The distribution of population-level genetic diversity of a species, i.e., its genetic population structure and the historical processes that formed that distribution (phylogeography), is a product of the species' life history and dispersal patterns, geographic history, climatic history, and chance (Avise 2007). An implicit assumption in the phylogeographic and population structure literature is that closely related species that occupy similar niches and occur in sympatry should maintain comparable population structure and share a similar population history. The comparative analysis of phylogeography and population structure of closely related, codistributed species provides a general framework within which this assumption can be tested. With carefully chosen model species, comparative phylogeography can highlight the relative contribution of life history vs species and population history factors and identify common processes that affect communities (Kelemen and Moritz 1999, Hodges et al. 2007). Therefore, comparative phylogeography of similarly distributed species can help reveal the existence of concordant or varying phylogeographic breaks between regions. The more taxonomically and ecologically diverse the selected species are, the more likely it is that common patterns in genetic structure will reflect general patterns for a wide range of faunas. Conversely, if codistributed species show different patterns of genetic structure, insights can be gained into species-specific characteristics, such as habitat specificity, dispersal ability, and resistance to barriers (Hodges et al. 2007). However, the effectiveness with which these results can be projected to other species of close evolutionary descent remains unknown. Closely related, preferably recently diverged, species should be examined to incorporate the role of the evolutionary template or constraints upon which intraspecific processes act in shaping the genetic structure of populations.
Comparative phylogeographic studies are still relatively rare and often compare very distantly related taxa (e.g., Taberlet et al. 1998, Zink et al. 2001, Carstens et al. 2005) to assess generality of patterns. However, several studies have shown that relatively closely related taxa have differing population histories (e.g., Hodges et al. 2007). Therefore, differences in patterns observed between distantly related taxa could be the result of their varied evolutionary history, as well as differences in ecology or distribution. Only a few studies have dealt with closely related taxa, including sibling ants (Pusch et al. 2006), bats (Berthier et al. 2006), scorpions (Salomone et al. 2007), mice (Michaux et al. 2005), and freshwater crayfish (Trontelj et al. 2005). These studies examined closely related, generally parapatric species pairs that occurred in sympatry only along secondary contact zones, and they often focused on patterns of introgression and hybridization. Michaux et al. (2005) examined 2 very closely related mice species that occurred in sympatry across a large part of their range in the Western Palearctic but that differed in their ecological plasticity.
Patterns of population structure in the European biota observed in comparative and single-species studies often are related to Pleistocene climate change. Effects of climate change on the biota have been shown for many plants and animals (see Hewitt [2004] and Schmitt [2007] for recent reviews). From recent studies, we know that aquatic organisms have a diverse array of population structures that suggest varied responses to past climate change (Wilcock et al. 2001, Damgaard 2005, Heilveil and Berlocher 2006, Pauls et al. 2006, Williams et al. 2006, Engelhardt et al. 2008, Lehrian et al. 2009). Often, these patterns differ from those observed in similarly distributed terrestrial taxa. For example, many aquatic species have occupied numerous microrefugia (e.g., Weiss et al. 2002, Pongratz et al. 2003, Trontelj et al. 2005, Verovnik et al. 2005, Pauls et al. 2006). The high numbers of refugia observed in the Iberian (e.g., Gomez and Lunt 2007) and Balkan Peninsulas (e.g., Krystufek et al. 2007) and in the Carpathian Mountains and basin (Botosaneanu 1973, 1975, Malicky 2006) indicate the importance of these regions as Pleistocene refugia and diversification centers.
We tested explicitly if the genetic population structure of closely related, partially sympatric species differs by comparing the population structure of Drusus discolor (Rambur 1942), Drusus romanicus romanicus Murgoci and Botosaneanu 1953, and Drusus romanicus ssp. meridionalisKumanski 1973. We used these species to test our hypothesis that closely related, sympatric, and ecologically similar species should show comparable patterns of population structure because they are putative sister species (Pauls 2004, Pauls et al. 2008a), occupy the same ecological niche (Botosaneanu 1959, Bohle 1983, Pauls et al. 2008a), and occur sympatrically in parts of the Carpathians and Balkan Highlands. For these species, we can assume that differences in population structure do not reflect a vastly different evolutionary history and do not result from occupation of different microhabitats. Rather, differences in genetic population structure could be caused, for example, by varying life histories or dispersal capacity.
Methods
Study system
Recent molecular phylogenetic studies have shown that D. discolor is the sister taxon to D. romanicus romanicus and D. romanicus meridionalis (Pauls 2004, Pauls et al. 2008a). The subspecies D. romanicus romanicus and D. romanicus meridionalis were described by Kumanski (1973) as geographically exclusive species with minimal but constant differences in the genital armature. We chose an inclusive definition of both D. discolor and D. romanicus because our main goal was to examine how 2 evolutionary lineages of Drusinae evolved in the region. Thus, we included all identified genetic lineages in both species. In D. romanicus, we also included both morphologically recognized subspecies because these subspecies might not reflect the whole extent of lineage diversity in the taxon, and similar lineages might have been found in D. discolor had comparable morphological studies been done. Thus, we combined D. r. romanicus + D. r. meridionalis in our genetic analyses and, hereafter, refer to them together as D. romanicus. We also used this opportunity to test the validity of the subspecies concept proposed by Kumanski (1973).
Drusus discolor is easy to distinguish based on adult male genital morphology from either subspecies of D. romanicus, whereas the females are difficult to distinguish based on morphology, and no distinctive morphological features are known to separate the larvae. Larvae of both species inhabit high-gradient, small- to mid-sized headwater streams and are generally restricted to altitudes >1000 m above sea level (asl) in the southern and eastern mountain ranges. The larvae of both species are filter-feeding carnivores and share unique adaptations in larval morphology and a unique feeding behavior (Botosaneanu 1959, Bohle 1983, Waringer et al. 2007a). Drusus discolor occurs in all European mountain ranges with elevations >800 m asl from the Iberian Peninsula in the west to the Rhodopi Mountains in the east, with exception of the Apennine Mountains in Italy. Drusus romanicus has a much more restricted range and has been recorded only from the Southern and Western Carpathians, the northwestern extensions of the Rhodopi Mountains (Fig. 1).
Specimen collections
We collected adults and larvae of D. discolor, D. r. romanicus, and D. r. meridionalis from sites across the Carpathians and the Bulgarian highlands in July and August 2003. We did additional sampling in the Romanian Carpathians in summer 2006 and spring 2007. We collected adults at lights or by sweeping riparian vegetation with a hand net, whereas larvae were picked individually from the substrate or captured with a dip net. We collected and analyzed 105 specimens of D. discolor, 61 D. r. romanicus, and 13 D. r. meridionalis from 22, 12, and 1 localities, respectively (Fig. 1, Table 1). We collected all specimens into 95% ethanol and stored them cool in the laboratory until deoxyribonucleic acid (DNA) was extracted. We deposited all voucher specimens in the collection of the Senckenberg Museum of Natural History in Frankfurt am Main, Germany.
Laboratory methods
We extracted whole genomic DNA using DNEasy Tissue Kits (Qiagen, Hilden, Germany) and amplified a 498-basepair (bp) region of mitochondrial DNA (mitochondrial cytochrome oxidase I [mtCOI]) via polymerase chain reaction (PCR) with the primers Jerry (Simon et al. 1994) and S20 (Pauls et al. 2006). Methods are described in detail in Pauls et al. (2006). We generated sequences with the PCR primers in the Field Museum Pritzker Laboratory for Molecular Systematics and Evolution, Chicago, Illinois (USA), and the Nano+Bio Zentrum at Kaiserslautern University, Kaiserslautern (Germany). We assembled and manually checked sequence traces with the program Seqman (DNASTAR, Madison, Wisconsin). We aligned sequences with Clustal W (Thompson et al. 1994) as implemented in Bioedit (Hall 1999).
Analyses
We followed 3 steps in our analytical procedure. First, we assigned larvae and females to either D. discolor or D. romanicus based on geographic distribution and haplotype comparisons with clearly morphologically distinct adult males. Second, we reconstructed the phylogenetic relationship between the mtCOI haplotypes of D. discolor and D. romanicus. Third, we used mtCOI haplotype divergence, frequency, and distribution to compare the genetic population structure of D. discolor and D. romanicus in the Carpathians and Bulgarian highlands, focusing especially on the regions of codistribution (Southern Carpathians and the northwestern Rhodopi Mountains).
Species assignment of larvae and females
We generated a neighbor-joining tree using the K2P model and estimated bootstrap support values for each node based on 1000 bootstrap replicates. This approach has been widely adopted for DNA barcoding studies (e.g., Ball et al. 2005, Zhou et al. 2007). The analyses were done in Paup* 4.0b10 (Swofford 2001). We followed the approach of Zhou et al. (2007) and considered larvae and females conspecific to an adult male if they either shared a haplotype with an adult male or they were nested within the species boundaries of adult males. If a larva or female was placed basal to the species boundaries based on adult males, the specimen was associated only if it was grouped in a highly supported, otherwise monospecific clade. We also used geography as a secondary criterion to verify associations, based on the fact that only 1 species has been recorded from many regions, e.g., only D. r. romanicus occurs in the Western Carpathians, and only D. discolor occurs in the Eastern and Northern Carpathians.
Phylogenetic relationship between D. discolor and D. romanicus
To examine the relationship between D. discolor and D. romanicus, we constructed a phylogenetic tree of all known mtCOI haplotypes, including those from other regions published in Pauls et al. (2006, 2008a) and 2 specimens each of Drusus chrysotus and D. muelleri as outgroup taxa. We did Bayesian phylogenetic analyses with the Markov chain Monte Carlo method (B/MCMC) in MrBayes 3.1 (Ronquist and Huelsenbeck 2003), assuming the GTR+I+G model of nucleotide substitution. We ran 2 parallel analyses with 12 chains each for 2 × 106 generations. We sampled trees every 100th generation. We discarded the first 1.5 × 106 generations as burn-in. We plotted the log-likelihood scores of sample points against generation time using TRACER 1.4 ( http://beast.bio.ed.ac.uk/Tracer) to ensure that stationarity was achieved after the first 1.5 × 106 generations by checking whether the log-likelihood values of the sample points reached a stable equilibrium plateau. From the remaining 10,000 trees, a majority-rule consensus tree with average branch lengths was calculated using the sumt option of MrBayes. Posterior probabilities were obtained for each clade.
We used the Shimodaira–Hasegawa (1999) (SH) test and expected likelihood weights test (ELW) (Strimmer and Rambaut 2002) to evaluate whether our data were sufficient to reject alternative topologies. Such topologies, which might not be significantly worse than the topology obtained, might be present in suboptimal trees not sampled or not present in the 50% majority-rule consensus tree of the MCMC sampling. We tested the following hypotheses: 1) monophyly of D. romanicus, 2) monophyly of D. discolor, 3) monophyly of D. r. romanicus. We did the SH and ELW tests with Garli (Zwickl 2006) and Tree-PUZZLE 5.2 (Schmidt et al. 2002). We did unconstrained and constrained maximum likelihood (ML) analyses in Garli using the GTR+I+G nucleotide substitution model. We compared a pair of trees, including the best tree agreeing with each of the null hypotheses, i.e., the constrained ML tree, and the unconstrained ML tree, in SH and ELW tests using the “user tree evaluation” option with accurate parameter estimation assuming the GTR model in Tree-PUZZLE 5.2 for each hypothesis.
Introgression and incomplete lineage sorting
Introgression (gene movement among species through hybridization; Avise 2004) and incomplete lineage sorting (an early stage of gene divergence, where genetic drift in one or both genetically distinct and isolated lineages has not progressed sufficiently to allow gene sorting) can confound relationships among species and render them not monophyletic. To examine the possibility of recent or ongoing gene flow between D. discolor and D. romanicus through introgression, we used the program IMa (Hey and Nielsen 2007). The isolation with migration model implemented in IMa has 6 demographic parameters that define the growth dynamics of a population. These parameters include 2 migration rates, one for each population. The IMa software estimates the posterior probability for each of the model parameters, fitting the Isolation with Migration model to the data. The program estimates parameters for a pair of closely related populations or species, so all sequences of each species were used in the analysis as a single population. We ran 3 replicate MCMC runs using the IMa software with different starting seed numbers to guarantee convergence of the sample. Based on exploratory runs using different settings and starting seeds, we set prior maximum values q = 50 for both species, migration in both directions to 2N based on FST (m = 0.66), and t to 50. With these settings, all model parameters returned unimodal distributions with clearly defined peaks.
Analysis of genetic population structure
Based on the clear differences in the genitalia of adult males and our specific assignment of larvae and females to those males, we divided the data into specific data sets for D. discolor and D. romanicus. We used these 2 data sets for all comparisons of genetic population structure. Because the number of individuals from single locations was sometimes very low, we pooled all samples from localities within a single mountain range under the assumption that dispersal within ranges was possible. Thus, our analyses were done at the scale of mountain ranges or regions.
We calculated MJ networks (Bandelt et al. 1999) for each individual species with the default setting in Network 4.1 (Fluxus Technology, Suffolk, UK). We subsequently color-coded the origin of each specimen carrying a given haplotype to illustrate haplotype distribution and to identify number of haplotypes and number of endemic haplotypes in each mountain range or region. We estimated gene diversity using the default settings in Arlequin 3.1 (Excoffier et al. 2005). We conducted analysis of molecular variance (AMOVA; Excoffier et al. 1992) to partition total variation among different geographic hierarchies. We used 2 different a priori grouping schemes for AMOVA: in grouping A, we analyzed sampling localities within and among mountain ranges; in grouping B, we analyzed mountain ranges within and among regions. This grouping followed major orogenic breaks in the region and the 500 m asl elevation, below which neither species has been found. The grouping also assumed that individuals could disperse among localities within mountain ranges. Preliminary tests of population differentiation between localities found no significant genetic population structure within mountain ranges (data not shown), supporting this assumption. To further investigate genetic population structure, we did Exact Tests of Population Differentiation (ETPD; Raymond and Rousset 1995) between mountain ranges and region. We did ETPDs in Arlequin with 15,000 steps in the Markov chain, of which we discarded 5000 as burn-in. We also did Mantel tests using 1000 permutations by correlating population pairwise FST values and pairwise geographic distances to examine whether an isolation-by-distance effect existed for both species. We calculated geographic distances with the Geographic Distance Matrix Generator (Ersts 2008) and the FST matrix and Mantel test in Arlequin 3.1 (Excoffier et al. 2005). We used tests of selective neutrality (Tajima's D: Tajima 1989a; Fu's F: Fu 1997) as a means to find evidence for potential recent increases in effective population size (Tajima 1989b).
Results
DNA sequencing resulted in 45 and 69 previously unpublished sequences for D. discolor and D. romanicus, respectively. With these sequences and previously published data from Graf et al. (2005) and Pauls et al. (2006, 2008a), we compiled a data set with 105 sequences for D. discolor and 74 sequences for D. romanicus. We found 23 haplotypes, 9 of them new, for D. discolor, and 12 haplotypes, 8 of them new, for D. romanicus. We deposited haplotype sequences in GenBank (see taxon labels in Fig. 2).
Phylogenetics and specimen assignment
Species assignment of larvae and females
Specimen assignment is summarized graphically in Fig. S1 (available online from: http://dx.doi.org/10.1899/08-100.1.s1 (10.1899_08-100.1.s1.pdf)). No haplotypes were shared between adult males of the 2 species, and the minimum number of bp differences was 8 (Bucegi Mountains). Based on our assignment criteria, 71, 53, and 10 larvae from the Eastern European mountain ranges were associated to D. discolor, D. romanicus romanicus, and D. r. meridionalis, respectively. Two females were associated to D. discolor, and 3 females were associated to D. r. romanicus. All but one of the associated specimens shared a haplotype with an adult male or were nested within a highly supported specific clade. The only exception was a larva (L1DDR162) that was directly basal to a male of D. r. romanicus and formed a distinct and highly supported clade (bs = 100) with other D. r. romanicus specimens. Furthermore, the larva was collected from the Apuseni Mountains, an area where only D. r. romanicus occurs. Thus, we associated this larva to D. r. romanicus.
Based on our final sampling, including associated specimens, we observed that both species occurred in sympatry in 3 localities: Piatra Alba (Bucegi Mountains), Rausor River (Retezat Mountains), and Banderishka River (Pirin Mountains). Both species occurred within the Fagaras Mountains, but not in the same sampling localities. With the exception of the Rausor River, the same was true in the Retezat Mountains. In the Northern and Eastern Carpathians, the Central Balkan Highlands, and the Rila Mountains, we found only D. discolor. In the Apuseni Mountains (Western Carpathians), we found only D. r. romanicus. These results of regional occurrence support distributions recorded in the literature (Kumanski 1988, Ciubuc 1993, Botosaneanu 1995). Both species are recorded from all mountain ranges of the Southern Carpathians, but we found them in the same streams in only 2 cases.
Phylogenetic relationship and potential migration between D. discolor and D. romanicus
A recent 3-gene molecular phylogeny of the subfamily Drusinae shows that D. discolor and D. romanicus are closely related sister taxa (Pauls et al. 2008a). The phylogenetic inference of all mtCOI haplotypes of D. discolor and D. romanicus from the current study and the adult genitalia of both species are shown in Fig. 2. The 2 species form a highly supported clade with respect to the selected outgroup species (clade A). The relationships within the discolor–romanicus clade remain largely unresolved. However, the 2 most basal, highly supported clades correspond to D. r. meridionalis (clade B) and D. r. romanicus (clade C). This result supports the subspecies concept proposed by Kumanski (1973). The only haplotype of D. r. romanicus that is not in either of these 2 clades is haplotype DR01, which was found in the Bucegi Mountain collections. It falls within D. discolor. Within D. discolor, there are 6 supported clades that correspond to geographic lineages as discussed in detail in Pauls et al. (2006). Our topology renders D. romanicus polyphyletic with 3 clades, one of which is nested within D. discolor. Drusus discolor is paraphyletic with respect to haplotype DR01. However, our 3 alternative hypotheses (monophyly of D. romanicus, D. r. romanicus, and D. discolor) could not be rejected in both SH and ELW tests (p > 0.05). The divergent D. romanicus haplotype (DR01) was carried by 2 males collected at the site. It diverges from all other D. romanicus haplotypes by 29 nucleotide changes (5.8% divergence), and from the closest D. discolor haplotype by 8 nucleotide changes (1.6% divergence). Thus, the 2 species can be differentiated clearly on the basis of both morphology and mtCOI haplotypes, but the relationship between D. discolor and D. romanicus cannot be resolved with the data at hand.
Introgression and incomplete lineage sorting
We used coalescent simulations to examine whether the isolation-with-migration model implemented in IMa could explain the lack of monophyly between D. discolor and D. romanicus. Maximum likelihood estimates of migration parameters revealed the highest likelihood at a migration rate of 0. Figure 3 shows the posterior distributions for migration rates and reveals the highest probability for no migration in either direction between D. discolor and D. romanicus.
Analysis of genetic population structure
Haplotype diversity, distribution, and divergence
Species-specific MJ networks are shown in Fig. 4. We identified 23 unique haplotypes in D. discolor (Fig. 4, Table 1). Forty-four steps were necessary to link all haplotypes into a single network. The MJ network revealed 3 diverged haplotype groups, separated by 13 and 14 nucleotide changes (2.61–2.81%). Group 1 (top left) consisted of 6 haplotypes found in the northwestern Rhodopi Mountains (Rila and Pirin Mountains) and Macedonia (Shar Mountains). Group 2 (top right) consisted of 3 haplotypes from the northwestern Rhodopi Mountains (Rila Mountains) and the Central Balkan Highlands (Stara Planina). The 14 remaining haplotypes formed group 3 (bottom). Haplotypes from this group were found in all regions except the Central Balkan Highlands and northwestern Rhodopi Mountains. Within group 3, 3 haplotypes (DD47, DD50, DD15) were shared among regions. The most common and central haplotype, DD47, was found in all regions of Carpathian Mountains and in Macedonia (Shar Mountains). Ten haplotypes were private to a single region. However, 8 of these were singletons, and it is unclear whether they are truly endemic to a certain region. In contrast, haplotype DD46 was found in 12 individuals from the Bucegi Mountains only and appears to be truly endemic. Categorization of haplotype endemism is summarized in Tables 1 and 2.
In D. romanicus, we identified 12 unique haplotypes, which were linked into a single network in 63 steps (Fig. 4, Table 1). Four geographically distinct lineages were observed. These differed from one another by 7 to 29 nucleotide changes (1.4–5.8%). Group 1 (top left) consisted of 3 haplotypes, and all individuals were from the northwestern Rhodopi Mountains (Pirin Mountains, white). Group 2 (bottom) contained all 3 individuals collected from the Bucegi Mountains in the southeastern Carpathians (light gray). These individuals all carried haplotype DR01. Groups 3 and 4 (top right) consisted of haplotypes found in the Western Carpathians (black) and the Southern Carpathians (light and dark gray), respectively. Group 4 was the only group with haplotypes that occurred in more than one region. Its central haplotype, DR02, was common in both the Fagaras Mountains (southeastern Carpathians, light gray) and across the southwestern Carpathians (dark gray). Thus, all but one haplotype (DR02) found in D. romanicus were private and appeared to be regionally endemic (Tables 1, 2).
Drusus discolor and D. romanicus showed similar levels of haplotype endemism at the mountain-range scale (74% and 75%, respectively). However, at the regional scale, all haplotypes except 1 (92%) were regionally endemic in D. romanicus, whereas 83% of haplotypes found in D. discolor were regionally endemic (Table 2).
General diversity indices are summarized in Table 3. In a direct comparison between regions where both species were sampled (southeastern and southwestern Carpathians, northwestern Rhodopi Mountains), D. discolor had more haplotypes than did D. romanicus in the southeastern Carpathians (5 vs 2) and northwestern Rhodopi Mountains (8 vs 3). Drusus romanicus had more haplotypes (4) in the southwestern Carpathians than did D. discolor (2). However, when corrected for number of specimens sampled, haplotype diversity between the 2 species was very similar in each of these regions. In contrast, gene diversity in D. discolor was higher than in D. romanicus in all 3 regions as a result of greater divergence between haplotypes in D. discolor and higher number of haplotypes within each of these regions.
Population differentiation and genetic population structure
We did the AMOVA analysis at 2 levels of geographic structuring (Table 4). All levels of testing showed highly significant levels of fixation (p < 0.002). In geographic grouping A, which partitioned variation among mountain ranges, the percentage of variation found at each level differed considerably for D. discolor and D. romanicus. In D. discolor, most of the variation was found between mountain ranges (73%), but considerable variation also was found within sites (∼19%), and some variation was found among localities within mountain ranges (∼8%). Almost all of the variation (∼97%) in D. romanicus was found between mountain ranges. Within locality (∼1.4%) and within mountain range (∼1.6%) variation contributed only minimally to the total variation. Grouping B examined variation among regions. Very similar values were found for all 3 hierarchical levels in D. discolor. In contrast, for D. romanicus, most of the variation was among regions (∼82%), but variation among localities within regions was large (∼17%). This shift in components of variation could have resulted from the high degree of differentiation and divergence between haplotypes found in the Fagaras Mountains and Bucegi Mountains within the southeastern Carpathians. In general, the AMOVA results reflected high levels of genetic population structure among geographically isolated mountain ranges and regions in the study area.
High degrees of population differentiation and genetic population structure in D. discolor and D. romanicus also were indicated by ETPD (Table 5). We used ETPD to examine genetic population structure among localities, mountain ranges, and regions. At all levels, the number of significant comparisons was almost 2× higher in D. romanicus as in D. discolor. In D. romanicus, 47% of comparisons were significant among localities, compared to 28% in D. discolor. Among mountain ranges, the number of significant comparisons increased for both species to 67% in D. romanicus and 36% in D. discolor. Among regions, all 6 comparisons (100%) were significant for D. romanicus, whereas only 13 of 28 comparisons (46%) were significant in D. discolor.
The Mantel test revealed a low but significant correlation between population pairwise FST and geographic distance for D. discolor (r = 0.319, p = 0.001). The correlation was higher for D. romanicus, indicating a slightly stronger isolation by distance effect (r = 0.522, p < 0.0001).
Population history
Both Tajima's D and Fu's F neutrality tests showed significantly negative departures from neutrality in the southwestern Carpathians for D. romanicus (D = −1.733, F = −3.324; Table 3). This case was the only one in both species at the regional and mountain-range scales for which both tests yielded significant negative results. Negative values also were found for D. discolor in the Northern Carpathians, but only Fu's F was significant (F = −1.797; Table 3). Neutrality tests were contradictory for D. discolor in Macedonia, showing significant negative values for Tajima's D (D = −1.467) and positive values for Fu's F (F = 5.304; Table 3).
Discussion
Assigning specimens
Use of DNA barcoding (Hebert et al. 2003) or DNA taxonomy (Vogler and Monaghan 2006) to associate different life stages and sexes is becoming more common and accepted. Aquatic ecologists and taxonomists are taking advantage of these advances to associate unknown larvae of aquatic insects to known species based on comparisons of larval and adult DNA (e.g., Willassen 2005, Graf et al. 2009a). In their work on Asian Hydropsychidae, Zhou et al. (2007) established a rigorous framework for testing species associations using 2 molecular loci. We associated larvae and females of D. discolor and D. romanicus to readily identifiable adult males with the association criteria set by Zhou et al. (2007), but based on a single gene. Several other investigators have used the section of mtCOI that we analyzed for D. discolor and D. romanicus to associate previously unknown larvae to male specimens in the subfamily Drusinae (Graf et al. 2005, 2009b, Waringer et al. 2007b, 2008). Thus, a possibility of misidentifications remains with our methods, given the limitations of a single marker, but we are confident that our associations are correct. We were able to associate life stages, even in a situation where sister taxa were not reciprocally monophyletic with respect to mtCOI, because all of our haplotypes were diverged by ≥1.6% from those of the other species.
Relationship of D. discolor and D. romanicus
Our study, previous molecular phylogenetic analyses (Pauls et al. 2008a), and the very close larval and female morphology clearly show that D. discolor and D. romanicus are sister taxa. However, the relationship of clades within these species cannot be resolved fully based on the mtCOI data at hand. Our hypothesis testing cannot rule out monophyly of both species, but none of the exploratory analyses (maximum parsimony, neighbor-joining, B/MCMC [results not shown]) ever showed monophyly of either species in topologies. Drusus romanicus haplotype DR01 was always nested within D. discolor. In general, absence of reciprocal monophyly does not necessarily warrant discarding currently recognized species (Carstens and Knowles 2007) because of commonly observed discordance between species trees and gene trees (Page and Charleston 1997). Lack of reciprocal monophyly can result from incomplete lineage sorting between recently diverged evolutionary lineages (e.g., Pollard et al. 2006), ongoing migration (introgression) between lineages, or historic introgression events that took place after the initial lineage split (e.g., Buckley et al. 2006, McGuire et al. 2007). Known tests for identifying incomplete lineage sorting use multilocus data and compare and contrast patterns of signal between markers based on the assumption that gene genealogies sort at different rates after a speciation event (e.g., Buckley et al. 2006, Linnen and Farrell 2007). Unfortunately, we cannot test for incomplete lineage sorting with the data at hand. However, in a case of incomplete lineage sorting, we would expect to find more mixing of lineages than was observed (e.g., Takahashi et al. 2001, Buckley et al. 2006, Carstens and Knowles 2007). The results from IMa suggest that no migration currently exists between the study species. Rather, the pattern for D. discolor and D. romanicus seems to be most consistent with a past introgression event, in which a rare interspecific mating took place between previously split D. discolor and D. romanicus at a site in the Bucegi Mountains. These lineages appear to have been diverging since the introgression event because the lineages we find in sympatry are diverged by 9 bp (1.8%). The interspecific mating that led to haplotype DR01 appears to have been a rare event because we found no evidence for introgression among any of the other sites where both species occur in sympatry. Unfortunately, we cannot fully resolve the situation with the single-locus data set.
Genetic population structure
Our study presents a direct test for a common genetic population structure and population history in closely related, partially sympatric, and ecologically similar species. The results show that both D. discolor and D. romanicus exhibit strong genetic population structure across the eastern European mountain ranges. Both species show similar patterns with: 1) many private or endemic haplotypes, 2) high levels of among-region and among-mountain range variation, 3) significant population differentiation among mountain ranges and regions, and 4) significant correlations between geographic distance and genetic differentiation.
However, the strength and detail of patterns differ between D. discolor and D. romanicus. The most striking differences are in gene diversity, regional differentiation, and lineage divergence. Within localities, gene diversity is higher in D. discolor than in D. romanicus because divergence between haplotypes is greater and the number of haplotypes within regions is higher in D. discolor, and almost no variability exists within populations of D. romanicus. Also, more haplotypes are shared among mountain ranges and regions in D. discolor than in D. romanicus, where only one haplotype is shared among regions. Together, these results suggest that more occasional gene flow, which helps maintain higher levels of genetic diversity, might occur among populations of D. discolor than of D. romanicus, and that isolation of lineages is older in D. romanicus than in D. discolor. In D. discolor, the 3 observed lineages are diverged by 2.6% to 2.8%. In contrast, the 4 lineages observed in D. romanicus are diverged by 1.4% to 5.8%. Also, in D. romanicus, each lineage is endemic to an individual mountain range. These results show that D. romanicus is more strongly regionally isolated and that lineages separated earlier than in D. discolor. This result is supported by minimal, but consistent, morphological differences between the Romanian and Bulgarian specimens of D. romanicus, which Kumanski (1973) recognized and formalized by describing D. r. meridionalis. However, another highly diverged (5.8%) lineage of D. romanicus, which is endemic to the Bucegi Mountains, exists and might represent a new subspecies. Only 2 adult males were collected from this lineage of D. r. romanicus. No morphological differences were detected between this and the other lineages, but the small sample size did not allow an exhaustive morphometric analysis.
Dispersal could play a role in shaping the differences in genetic population structure of D. discolor and D. romanicus. Based on the lower levels of divergence among lineages, shared haplotypes, and slightly lower rates of genetic differentiation, we might expect D. discolor to be the better disperser of the 2 species. For example, the much larger range across Europe occupied by D. discolor could be the consequence of better dispersal. If D. discolor is truly more vagile than D. romanicus, we would expect evidence for recent range expansion. However, our neutrality test results showed such a pattern only for D. romanicus in the southwestern Carpathians, where a common haplotype is shared between all southern Carpathian ranges except for the Bucegi Mountains. Thus, no conclusive evidence shows that D. discolor is a better disperser. Alternatively, the larger range of D. discolor could also be the result of competition between the 2 species. Both species apparently occupy the same niche and occur in partial sympatry; thus, D. discolor might have outcompeted D. romanicus, enabling it to occupy a much larger range.
The comparative aspect
Our study shows that 2 very closely related, partially sympatric caddisfly species with very similar ecologies differ in their genetic population structure. Thus, neither a close relationship between taxa nor similar ecologies and distributions appears to be a good predictor for a common population history. Comparative phylogeography studies are becoming more common (Hodges et al. 2007, Joger et al. 2007), but most studies comparing genetic population structure either examine more distantly related taxa with similar ecologies or distributions (Joger et al. 2007, Lehrian et al. 2009) or study sibling taxa with largely allopatric distributions (Berthier et al. 2006, Pusch et al. 2006, Salomone et al. 2007). Thus, these studies ignore the effects that a different evolutionary ancestry or local differences in Pleistocene climate change might have had on the examined population histories. For example, Lehrian et al. (2009) showed that 2 distantly related caddisfly species with a similar distribution range had very different patterns of genetic population structure in Central Europe. The authors interpreted these differences to mean that one species (D. discolor) survived Pleistocene glaciations in several refugia in the periglacial mountains of Central Europe, whereas the other species (Hydropsyche tenuis) recolonized Central Europe, presumably postglacially, from a single or few southern European refugia. However, they were unable to rule out the possibility that the differences also could be attributed to different evolutionary histories of the 2 different families of caddisflies they examined. Salomone et al. (2007) examined 3 closely related species of scorpions in Italy. One species occurs in southern Italy and Tuscany, one in central and northern Italy, and the third inhabits central and northern Italy and Elba. Each of these regions was affected differently by Pleistocene climate change (Salomone et al. 2007), making it difficult to discern whether the species would have reacted similarly if subjected to the same climatic shifts. Pusch et al. (2006) examined the genetic population structure of 2 parapatric sibling ant species in a hybrid zone in Central Europe. Both species have similar ecologies, were exposed to similar Pleistocene climatic shifts, and have low levels of mtDNA sequence and allozyme signal variation. This pattern suggests that both species experienced a major bottleneck caused by very similar shifts in Pleistocene climates. Trontelj et al. (2005) showed that 2 sister taxa of freshwater crayfish in the genus Austropotamobius had similar patterns of genetic population structure and lineage divergence. In part as a result of human translocations, Austropotamobius pallipes is widespread across southwestern and Central Europe and also inhabits the region southeast of the Alps and the coastal Dinaric mountains. Austropotamobius trifidus, on the other hand, occurs in the upper Rhine Valley, in the southeastern Alpine region, along the Dinaric coastal mountains, and in the Rhodopi Mountains. Thus, both species occur in largely overlapping distribution areas southeast of the Alps and in the coastal Dinaric mountains. In contrast to the situation with D. discolor and D. romanicus, no large differences exist between these species in nucleotide diversity or average nucleotide divergence of mtCOI haplotypes. Both species harbor sufficiently divergent haplotype groups to indicate multiple refugia (2 in A. pallipes, 3 in A. trifidus) on the Balkan Peninsula. Many highly diverged haplotypes occur in A. trifidus, which occurs in the southern Balkan Peninsula. This pattern is unlike the patterns found in either species in any other region, but it is similar to our observations for D. discolor. Michaux et al. (2005) also found strikingly different patterns of genetic population structure in 2 closely related mice species that are sympatric across a wide range in Europe. These examples and our own study show that a similar evolutionary history, ecology, and distribution can sometimes lead to similar patterns of population history, but this outcome is not always the case and should not necessarily be expected.
The Carpathian and Balkan ranges—centers for diversification and Pleistocene refugia
The importance of the Carpathians and Balkan Highlands as centers of diversification and endemism has been discussed for many groups of organisms, including plants (e.g., Trigas et al. 2007, Hajek et al. 2008), mammals (e.g., Krystufek and Griffiths 2002, Kotlik et al. 2006), reptiles (Bohme et al. 2007), leeches (Sket and Trontelj 2008), spiders (Deltschev 1999), and a wide array of insects including butterflies (Schmitt and Haubrich 2008, Varga 2008) and caddisflies (e.g., Botosaneanu 1973, 1975, Mey and Botosaneanu 1985, Pauls et al. 2006, Balint et al. 2008, Previšic et al. 2009).
According to Botosaneanu (1975), 15% of the Carpathian caddisfly fauna is endemic. Particular centers for caddisfly diversity hosting several endemic species in the Carpathians are the Tatra Mountains (e.g., Drusus doehleri, Allogamus starmachi), the Southern Carpathians (Retezat, Fagaras, and Bucegi Mountains; e.g., Rhyacophila fagarasensis, Rhyacophila cibinensis, Allogamus dacicus), and the Apuseni Mountains (e.g., Drusus buscatensis, Rhyacophila orghidani). The Balkan mountain ranges, in particular the Rhodopean zone, have a large number of endemic caddisfly species, which are predominantly cold-adapted Limnephilidae and Rhyacophilidae that inhabit montane (>700–800 m asl) streams and lakes (Kumanski 2005). Within the genus Drusus, 31 species and subspecies are endemic to the Carpathians or Balkan Peninsula (Botosaneanu 1975, Pauls 2004, Kumanski 2005); this pattern demonstrates the importance of this region of the western Palearctic as a center of origination, particularly for caddisflies (Kumanski 2005).
Although the region is known for its high diversity and endemism, the results of our study warrant a discussion of the Bucegi Mountains, a mountain range of particular importance that is not often highlighted. Despite the geographic proximity of the Bucegi Mountains to other mountains sampled in our study (e.g., Fagaras Mountains), clear genetic isolation was found in both D. discolor and D. romanicus between specimens from the Bucegi Mountains and specimens from other regions in the Carpathians. The Bucegi Mountains are in the southeastern section of the Carpathian range and have peaks that are ≥2500 m asl in altitude. Because of their central position within the Carpathians, the Bucegi Mountains support a large number of Carpathian biotic elements. Their topography is diverse, including high-altitude glacial valleys, steep slopes to the north, and numerous caves, gorges, and sinkholes that are indicative of the calcareous geology. Together, these features provide good preconditions for diversification through allopatric speciation and have allowed many microendemic taxa to evolve. These endemics span many organism groups and include lichens (e.g., Polyblastia butschetschensis, Lecanora verrucosa var. bucegica, and Bucegia romanica; Mohan et al. 1993), higher plants (e.g., Thesium kernerianum, Hesperis oblongiflora, Saxifraga demissa, Astragallus australis ssp. bucegi, Poa molineri ssp. gelida, Festuca bucegiensis; Mohan et al. 1993), amphipods (e.g., Niphargus carpathicus cavernicolus; Bleahu et al. 1976), and stoneflies (e.g., Isoperla carpathica; Kis 1971).
The Balkan Peninsula, and more recently the Carpathian Mountains and basins, were identified as major refugial areas throughout Pleistocene cold phases (recently reviewed by Hewitt 2000, 2004, Schmitt 2007, Varga 2008). Some of the diversification within the region has been attributed to Pleistocene lineage divergence in independent localities within the Carpathian range (Botosaneanu 1973, Mey and Botosaneanu 1985, Stauffer et al. 1999, Schmitt et al. 2007) or between Balkan mountain ranges (Krystufek et al. 2007). More often, there is indication for preglacial lineage divergence (e.g., Hofman et al. 2007) and that lineages remained independent throughout much of the Pleistocene despite recurring population regression and expansion (Pauls et al. 2006).
It is difficult to draw phylogeographic conclusions for D. discolor and D. romanicus in the studied region with the data from our study. For example, the observed lineages are too diverged in both species to allow us to conduct phylogeographic inference using nested clade analysis (Templeton 1998) to explore potential hypotheses of population history. Also, with only a single marker, the data allow little room for coalescent simulations of population histories, which would allow explicit testing of phylogeographic hypotheses. Despite these drawbacks, patterns similar to those observed in D. discolor and D. romanicus—genetic differentiation and divergence among the Northern, Eastern, Southern, and Western Carpathians and the Balkan Peninsula—also have been observed in plants (e.g., Mraz et al. 2007, Ronikier et al. 2008), other caddisflies (Botosaneanu 1973, Mey and Botosaneanu 1985, Balint et al. 2008), and other insects (e.g., Schmitt et al. 2007), and allow a cautious comparative interpretation of our data. In all of these studies, the patterns were interpreted to indicate several independent Pleistocene refugia within the Carpathians or Balkan Highlands. Pauls et al. (2006, 2008b) reported that at least 2 independent refugia for D. discolor existed in the Balkan Mountains. This conclusion is supported by the 2 diverged lineages of D. discolor found in both the previous and current studies. Clear genetic population structure also exists within the Carpathians, but lineage divergence is much lower. The range-wide sampling of D. discolor indicates that the Carpathians served as a refugia for D. discolor (Pauls et al. 2006, 2008b). However, no evidence exists for recent increases in effective population size, so it is unclear whether extant genetic population structure within the Carpathians is the result of lineage divergence in independent Carpathian glacial refugia and postglacial secondary contact or of recent diversification. In D. romanicus, lineages clearly have diverged among the Southern Carpathians, the Western Carpathians, and the Bucegi Mountains. The lineage in the Balkan Highlands is also highly divergent. Together, these data suggest that D. romanicus could have persisted in 4 independent refugia, 3 in the Carpathians and 1 in the Balkan Highlands, for several glacial cycles.
Acknowledgments
We thank Krassimir Kumanski and Stoyan Beshkov (Sofia) for collecting much of the material used in our study with SUP in 2003. We also thank all the other colleagues who provided us with specimens, without which this study would not have been possible (Table 1). We thank Peter Tiffin (St Paul) for assisting us with the interpreting the IMa results, and Sharon Jansa and Keith Barker (St Paul) for comments on the study. We appreciate the helpful comments from Ralph Holzenthal (St Paul), Karl Kjer (New Brunswick), Pamela Silver (Erie), and an anonymous referee on an earlier version of this manuscript. Part of this work was done with the resources of the Computational Biology Service Unit from Cornell University, which is partially funded by Microsoft Corporation. Our study was financially supported by 2 German Academic Exchange Service (DAAD) Awards to SUP, and the research funding programme “LOEWE—Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz” of Hesse's Ministry of Higher Education, Research, and the Arts. This study is part of the outcome of the German Science Foundation research grant (DFG-Ha3431/2-1) awarded to PH and SUP.
Literature Cited
Figure 1.
Sampling locations of Drusus discolor (Rambur 1942) and Drusus romanicus Murgoci and Botosaneanu 1953. Shading corresponds to that used in Fig. 4. Mts = mountains.
Figure 2.
Bayesian Markov chain Monte Carlo inference for 81 mitochondrial cytochrome oxidase I (mtCOI) haplotypes from Drusus discolor (Rambur 1842) (HDd; in black), Drusus romanicus romanicus Murgoci and Botosaneanu 1953, and Drusus romanicus meridionalisKumanski 1973 (HDr; both in gray). The 50% majority-rule consensus tree from 10,002 trees sampled in 2 parallel runs with 5 heated and 1 cold chain each is shown. Tree space was searched for 5,000,000 generations. Bold branches indicate posterior probabilities ≥0.95. The tree is rooted with mtCOI haplotypes from the 2 putative closest relatives, Drusus chrysotus Rambur 1842 and Drusus muelleri MacLachlan 1868 (Pauls et al. 2008a). * indicates haplotypes that are present in the study region and were used in regional analyses; ˆ indicates new haplotypes found in this study. Taxon labels indicate species-specific haplotype numbers used in Fig. 4 and in Pauls et al. (2006). Taxon labels are followed by GenBank accession number for each haplotype. Labeling of clades (A, B, C) refers to clade descriptions in the text. The inset shows genitalia of adult male specimens of D. discolor and D. romanicus romanicus.
Figure 3.
Marginal posterior probability distribution for migration rate estimates under the “Isolation with Migration” model from Drusus discolor toward Drusus romanicus and D. romanicus toward D. discolor.
Figure 4.
Median-joining networks of haplotypes from Drusus discolor (Rambur 1942) (A) and Drusus romanicus Murgoci and Botosaneanu 1953 (B). Shading corresponds to origin of specimens carrying a given haplotype. Circle diameter is relative to number of specimens carrying a haplotype. Hatch marks indicate number of basepair changes between haplotypes. Shading corresponds to that used in Fig. 1. Mts = mountains.
Table 1.
Sampling localities and number of individuals of D. discolor used in this study. The specimens are organized in a geographic hierarchy reflecting the different hierarchical groupings we examined in the analysis of molecular variance (AMOVA). NR = number of individuals per region, NM = number of individuals per mountain range, NP = number of individuals per sampling locality. H = haplotype, R = regionally endemic haplotypes, M = mountain range endemic haplotypes. Country codes: BG = Bulgaria, MZ = Macedonia, PL = Poland, RO = Romania, SK = Slovakia, UA = Ukraine.
Table 1.
Continued.
Table 2.
Assessment of haplotype endemism for Drusus discolor and Drusus romanicus.
Table 3.
Gene and nucleotide diversity. * = significant Tajima's D (p < 0.05), ** = significant Fu's F (p > 0.02) after Excoffier et al. (2005).
Table 4.
Results of the analysis of molecular variance (AMOVA; Excoffier et al. 1992) for Drusus discolor and Drusus romanicus.
Table 5.
Results of Exact Tests of Population Differentiation for Drusus discolor and Drusus romanicus.
[1] This study is dedicated to the memory of Krassimir Kumanski who passed away in July 2006, 3 y after collecting much of the material used in this study with SUP.
[2] E-mail addresses: pauls497@umn.edu
