Rapid ecological isolation and intermediate genetic divergence in lacustrine cyclic parthenogens
© Costanzo and Taylor. 2010
Received: 30 October 2009
Accepted: 5 June 2010
Published: 5 June 2010
Skip to main content
© Costanzo and Taylor. 2010
Received: 30 October 2009
Accepted: 5 June 2010
Published: 5 June 2010
Ecological shifts can promote rapid divergence and speciation. However, the role of ecological speciation in animals that reproduce predominantly asexually with periodic sex and strong dispersal, such as lacustrine cladocerans, is poorly understood. These life history traits may slow or prevent ecological lineage formation among populations. Proponents of the postglacial ecological isolation hypothesis for Daphnia suggest that some species have formed postglacially in adjacent, but ecologically different habitats. We tested this hypothesis with ecological, morphological, and multilocus coalescence analyses in the putative lacustrine sister species, Daphnia parvula and Daphnia retrocurva.
Daphnia parvula and D. retrocurva showed strong habitat separation with rare co-occurrence. Lakes inhabited by D. parvula were smaller in size and contained lower densities of invertebrate predators compared to lakes containing D. retrocurva. In the laboratory, D. retrocurva was less vulnerable to invertebrate predation, whereas D. parvula was less vulnerable to vertebrate predation and was smaller and more transparent than D. retrocurva. The species are significantly differentiated at mitochondrial and nuclear loci and form an intermediate genetic divergence pattern between panmixia and reciprocal monophyly. Coalescence and population genetic modelling indicate a Late or Post Glacial time of divergence with a demographic expansion.
Despite their young age and mixed breeding system, D. parvula and D. retrocurva exhibit significant ecological and genetic divergence that is coincident with the formation of deep temperate glacial lakes. We propose that predation may have facilitated the rapid divergence between D. parvula and D. retrocurva and that intermediate divergence of aquatic cyclic parthenogens is likely more common than previously thought.
Divergent selection for different habitats has long been proposed as a contributing factor to the evolution of species [1–4]. Regardless of geographic associations, ecological speciation from strong single selective agents, such as predation, or from the additive effects of multiple selective agents (i.e., the niche dimensionality hypothesis) now have empirical support [5–10]. A remaining theoretical gap is a detailed understanding of the role of the breeding system in speciation [11–14]. Although many eukaryotic taxa reproduce both sexually and asexually , there are still few empirical studies of speciation beyond strictly sexual taxa [14–17]. Because many organisms with mixed breeding systems are excellent dispersers [18–21], breeding systems with even a small amount of sexual reproduction should suffer homogenization of ecologically differentiated lineages [11, 12, 15, 22]. Thus, taxa with mixed breeding systems coupled with high dispersal should possess a weak capacity for rapid ecological radiation compared to strict asexuals and sexuals [11, 23–26]. Others, however, have proposed that priority effects where early colonization offers an advantage and divergent selection are frequently strong enough to overcome the homogenizing effects of such mixed breeding systems [22–27].
Cladocerans contain several candidate groups for detailed empirical study of ecological and genetic divergence in a system characterized by both sexual and asexual reproduction (cyclic parthenogenesis). These microcrustaceans often inhabit insular freshwater systems (i.e. lakes, ponds, reservoirs) with a mosaic of selective regimes through many environmental factors such as ephemerality, temperature, nutrient availability, competition, and predation [28–30]. For example, a common contrast occurs between lakes with visual vertebrate predators that can select for smaller prey species, and lakes dominated by gape-limited invertebrate predators that favor cladoceran communities dominated by large or helmeted species [31–36], as prey vulnerability is a function of prey size relative to predator size. Predation has been proposed to drive evolutionary divergence in morphology and behaviour [32, 37–41] and speciation in cladocerans [23, 25, 26, 41–43]. Still, there is little knowledge of ecological speciation among sister species of lacustrine cladocerans, where sexual reproduction is relatively infrequent.
Daphnia parvula lacks the anterior retrocurved helmet (Fig. 1) that appears to be controlled by polyploid cells in the head region of D. retrocurva . The helmet undergoes cyclomorphosis in D. retrocurva [32, 56–60], which can be enhanced by cues from invertebrate predators [32, 58–60] and laboratory studies also show decreased vulnerabilities of D. retrocurva to invertebrate predation with increased helmet size . Importantly, multigenerational lab culture experiments failed to convert the head shape of D. retrocurva into the non-helmeted head shape of D. parvula [56, 62, 63]. Although the body size and tail spine length in D. parvula can be influenced by predator cues, attempts to induce helmet production have also been unsuccessful , suggesting the morphological differences between these two species likely has a genetic component. Additionally, although D. retrocurva is a classic example of seasonal polymorphism, experiments reveal that phenotypic plasticity also fails to explain the morphological differences between the two species [62, 63].
Although there are no studies that specifically address ecological differentiation between D. retrocurva and D. parvula, there is substantial evidence that the interaction of predation with helmet phenotypes is evolutionarily and ecologically relevant. Kerfoot and Weider  for example, show that the abundance of D. retrocurva and larger size classes of helmets in D. retrocurva, are positively associated with the increased abundance and size of the invertebrate cladoceran predator, Leptodora, and reduced levels of planktivorous fish over time. Likewise, following the introduction of Leptodora in a lake, a shift from D. parvula (without the helmet) was observed to helmeted or large Daphnia , suggesting a role for predation as a selective agent for helmets. In contrast, a sixteen-year study  found repeated smooth, continuous transitions between large bodied Daphnia and small Daphnia (including D. parvula) associated with planktivorous fish density. Additionally, the abundance of D. retrocurva in glacial lakes is positively correlated with the abundance of invertebrate predators, and the abundance of introduced D. parvula is negatively associated with rates of invertebrate predation during a season . Taken together, the translocation and natural experiments are consistent with divergent selection from predation playing an isolating role, which could also affect genetic differentiation. However, Lukaszewski et al.  found evidence that ecological conditions beyond predation can prevent the successful colonization of D. retrocurva into a lake containing D. parvula.
Isolation affects population genetic patterns in a time dependent fashion. Omland et al.  defined "intermediate divergence" as the stages between panmixis and reciprocal monophyly. The lineage divergence continuum, as visualized by networks, begins with gene frequencies diverging with only older internal haplotypes shared (less frequent haplotypes arising since isolation are private), followed by a lack of haplotype sharing. Eventually, population-specific clades (monophyly or paraphyly) are formed. In post-glacial isolation, which has been proposed by Brooks for D. parvula and D. retrocurva , such intermediate divergence is often produced and haplotype sharing is limited to the older central haplotypes of the networks [67–69]. In contrast, if isolation is caused by multiple glacial refugia during the late Pleistocene, patterns of multiple geographic subclades along with multiple demographic expansions is predicted as observed in other Daphnia [67, 68].
Here, we specifically test for ecological and genetic divergence between D. parvula and D. retrocurva. Ecological speciation predicts that populations of D. parvula and D. retrocurva whose ranges overlap will be associated with ecologically different habitats. If predation is a major selective agent, then we should detect differences in the relative vulnerabilities to habitat-specific predators. If the morphotypes of "retrocurva" and "parvula" are isolated, then we expect to detect significant genetic differentiation in the zone of overlap. Further, under the postglacial isolation hypothesis, we expect intermediate divergence (lack of monophyly), a lack of multiple demographic expansions, and divergence times from coalescence modelling to be Late or Post-Glacial. We test these objectives through data culled from pre-existing databases, laboratory predation experiments, and population genetic analyses.
Discriminant Analyses on abiotic variables.
Canonical variate correlation coefficients
Calculated Alkalinity (ueq/L)
Chlorophyll A (ug/L)
Ionic Strength (M)
Total Nitrogen (ug/L)
Total Phosphorus (ug/L)
Mean Secchi dish depth (m)
Total Suspended Solids (ug/L)
Mean Lake Depth (m)
Lake Volume (m 3 )
Lake Area (ha)
% Watershed in human disturbed land
Wilks' Lambda P
Of the 97 lakes analyzed from the Department of Natural Resources (DNR) in Wisconsin, there were 16 lakes with only D. parvula, 79 lakes with only D. retrocurva, and two lakes with both species. In both co-occurrence lakes, one species dominated while the other was rare or only present in a single sampling period throughout the year.
Discriminant Analyses on invertebrate predator abundances.
Canonical variate correlation coefficients
Wilks' Lambda P
Prey vulnerabilities to tested invertebrate and vertebrate predation.
Invertebrate Predation Experiment
Vertebrate Predation Experiment
In the vertebrate predation experiment, there were significant effects of treatment (control vs. predator) and treatment by species interaction for the proportion on prey missing (Table 3). There was no significant difference between the proportion of prey missing for D. parvula and D. retrocurva in the control treatments, but D. retrocurva suffered higher mortality by emerald shiners in the predator treatments than did D. parvula (Fig. 4b). The predation rate of the fish on D. retrocurva was also higher (1.085 ± 0.0296) than that on D. parvula (0.638 ± 0.033).
Detection of exclusive ancestry.
(HSP90, F6F12, G6G12)
Source of Variation
All 4 markers
Among populations within species
All 4 markers
ND2 HSP90 All 4 markers
All 4 markers
Our findings illustrate rare co-occurrence of D. parvula and D. retrocurva in the field with strong spatial segregation in ecologically differentiated lakes supporting ecological isolation. Priority effects and dispersal limitation can contribute to habitat isolation [71, 72], but do not predict a pattern of distribution associated with lake type. Moreover, both species appear to be strong dispersers [47, 64, 73] and unidirectional translocation studies have failed to result in colonization . The rare co-occurrence of parvula-like morphotypes with retrocurva described by Brooks  remains mysterious. We failed to detect parvula-like specimens in multiple samples from the same embayment of Lake Ontario examined by Brooks. The rare parvula-like morphs in glacial lakes could be the result of recent inflow from shallow adjacent waters or transient ex-ephippial clones of parvula.
What are the ecological traits associated with isolation? Kerfoot and Weider  have elegantly shown that evolution of the helmet size (defensive structures that are the main morphological difference between D. retrocurva and D. parvula) can be attributed to historical changes in invertebrate predation. Here we show that this process may also contribute to the spatial segregation of these two species in nature, as D. retrocurva inhabit larger lakes that contain higher abundances of invertebrate predators relative to lakes where D. parvula is found. Although the data collected did not permit us to examine if fish, or vertebrate predators discriminated between the habitats encountered by these two species, the abundance of D. parvula has been found to be positively associated with the abundance of planktivorous fish . We also were unable to make a direct association between the abiotic and biotic factors compared in this study since these data were obtained from different sources, yet deeper aquatic habitats can have higher abundances of invertebrate predators of Daphnia because depth affords a refuge from visual predators for conspicuous invertebrate predators . Unexpectedly, the lake data failed to reveal an association between the abundance of Leptodora and Daphnia species. One possible explanation for the lack of association is the difficulty in quantifying abundances of vertically migrating Leptodora in daylight hauls [75–77], yet the abundances of D. retrocurva and Leptodora have been found to be positively associated in other studies .
We also found that D. parvula has an advantage relative to D. retrocurva in the presence of visual, vertebrate predation which may be a consequence of the smaller body size and greater transparency of D. parvula (Fig. 5). Additionally, D. retrocurva has an advantage over D. parvula in the lab to the invertebrate predator tested. Thus, the present study and existing translocation and longitudinal evidence suggest that D. parvula and D. retrocurva have a selective advantage to different habitat-associated predators. It is possible that predation may promote "immigrant inviability" between these two species, or selection against migrants between locally adapted populations, which is a direct result of natural selection reducing gene flow [78, 79]. The differing predation regimes likely contribute to the spatial segregation of these sister species, but additional local processes may to contribute as well in this system . Our study found that abiotic chemical variables (i.e. nutrients, conductivity) did not discriminate between the two species' habitats, yet other processes that were not included in this study could also contribute to their habitat segregation such as competition, the interactive effects of multiple predators, priority effects, and dispersal [80–83]. Further studies are needed to evaluate the role these other factors may play in facilitating habitat isolation between D. parvula and D. retrocurva, yet evidence suggests that predation is a key factor in determining zooplankton community composition and can reduce the importance of other factors influencing the establishment of Daphnia populations [81, 84].
Our results support postglacial genetic divergence between D. parvula and D. retrocurva despite a mixed breeding system and high dispersal. The networks, single demographic expansion patterns, and divergence estimates based on nuclear and mitochondrial DNA are consistent with intermediate divergence and postglacial isolation. This postglacial timing of divergence is consistent with the origins of the deep temperate glacial lakes [85, 86] that constitute the main habitat of D. retrocurva. We note that the findings of intermediate divergence do not rule out actual or potential gene flow between D. retrocurva and D. parvula. The AMOVA reveals a relatively low proportion of the total genetic variation is found between the two species and genetic divergence associated with ecological processes is often characterized by intermediate divergence with ongoing gene flow .
Although Daphnia are known to have high dispersal and colonization rates , the fixation indices suggest strong genetic differentiation among intraspecific populations (Table 5), suggesting natural selection may also be operating within species as well. Therefore, selection against migrants may lead to a reduction in gene flow in this system on two scales, between two species and the two larger classes of habitat types they reside in, and within a species and the individual local habitat each population resides in. According to island population models, as long as the island population size increases more rapidly than the number of migrants, the per-generation genetic contribution of migrants should decrease and an equilibrium is reached  which is likely in Daphnia. Their cyclic parthenogenic reproduction, rapid population growth, and priority effects supported by a large resting bank coupled with local adaptations may allow residents to maintain large populations relative to incoming migrants. Many studies illustrate strong, local adaptations in Daphnia populations and the reduced success of invaders [21, 72, 81], this high genetic differentiation among intraspecific Daphnia populations often leads to a monopolization effect [67, 72, 88, 89] in which predation, competition, reproduction, ontogeny, and tolerance to eutrophication have been implicated to promote [72, 88–90].
Conflicting models have been proposed whereby mixed breeding systems either slow  or accelerate divergence [18, 23, 26]. In our study, despite their young age, D. retrocurva and D. parvula illustrate ecological divergence with respect to their habitats, morphology, and predator vulnerabilities suggesting possible selective advantages to each in their respective environments. Both species typically have a mixed breeding system with relatively infrequent sex [42, 44, 45, 47] and illustrate strong dispersal capabilities [47, 64, 73]. Although some believe that populations with mixed breeding systems are less likely to yield discrete or distinct groups relative to populations with strictly sexual or asexual reproduction , this study illustrates two groups that have diverged into two distinct lineages and maintained this divergence with a mixed breeding system coupled with high dispersal, and more importantly over a relatively short period of time. Lynch and Gabriel  propose that a single generation of sex can reveal 50-75% of the genetic variance hidden by asexual reproduction. Therefore, reproduction marked by cyclic parthenogenesis could allow such rapid divergence through the combined benefits of rapid population growth rates that asexual reproduction offers along with the release of genetic variation in the few episodes of sexual reproduction . We postulate this mixed breeding system permits such rapid postglacial divergence of D. parvula and D. retrocurva, with selection against migrants (immigrant inviability) maintaining this divergence by reducing the homogenizing effects of dispersal.
We present evidence for rapid ecological, morphological, and genetic divergence of sister species of lacustrine cyclic parthenogens. The species differ in genetically-based defensive structures and predation regimes. Predation has been implicated as a selective force promoting divergence among taxa in only a few animal lineages. The present study suggests a role for predation in facilitating habitat isolation in young sister species with a mixed breeding system. However, multifarious selection can bring about reproductive isolation more effectively than selection from a single factor, and the habitats of D. retrocurva and D. parvula differ in numerous ways that could affect selection. Although predation has been implicated as a key factor in determining community composition and establishment of Daphnia [81, 84], other factors need to be further evaluated. Our results do indicate that a better understanding of the evolution of the numerous presumed intraspecific defensive-morphotypes of cladocerans is needed. Intermediate divergence of aquatic cyclic parthenogens is likely more common than previously thought.
We obtained occurrence data from the United States Environmental Protection Agency's (EPA) Environmental Monitoring and Association Proposal (EMAP), http://www.epa.gov/emap/. These data were collected from 64 lakes and ponds in the northeastern United States in the months of July through September from 1991 to 1994. From this study, data were extracted on the presence or absence of D. parvula or D. retrocurva within each habitat, along with physical and chemical measurements of: ionic strength (M), pH, total nitrogen (μg/L), total phosphorus (μg/L), trichomatic chlorophyll A (μg/L), total suspended solids (μg/L), secchi depth (m), calculated alkalinity (meq/L), turbidity (NTU), average depth of lake (m), lake volume (m3), lake surface area (ha), and % of watershed in human disturbed land. The methods used to measure each variable are described at http://www.epa.gov/emap/. We coded the presence/absence data in each site for both D. parvula and D. retrocurva. There was one lake in the EPA abiotic dataset that included both species, and in this lake D. parvula was rare. Data from this co-occurrence lake were removed from the analysis (see below) due to the small sample size.
Separate data collected by the Wisconsin Department of Natural Resources (DNR)  from 1973-1974 from 99 lakes throughout Wisconsin were used to test for differences in invertebrate predator assemblages across the species' habitats. We were unable to assess differences between vertebrate predator assemblages (fish) due to incomplete or missing data. Each lake was sampled four times that year (between the months of April through November) where zooplankton and invertebrate predator abundances were recorded. These data include ranked abundances of the following invertebrate predators: the copepods Acanthocyclops vernalis, Epischura sp., and Mesocyclops edax; the larval dipterans Chaoborus flavicans, C. punctipennis, and Chaoborus sp.; and the cladoceran Leptodora kindti. The abundance data were ranked according to the absolute abundances found in the field as follows: 0 = absent, 1 = rare (3 or less), 2 = uncommon (4-50), 3 = common (51-1,000), and 4 = abundant (> 1000). We coded the presence/absence of D. parvula and D. retrocurva. Since the timing of the four periods sampled varied across lakes, the mean values for each lake were used in the analyses. This dataset had two lakes where the species co-occurred (never during the same time period), which were removed from the analysis (see below) due to the small sample size.
Because abiotic and biotic variables are often correlated , the datasets were analyzed with a Discriminant Analysis (DA) in SPSS (ver. 11.5, 2006). Two separate DA's were run, one including the abiotic data (EPA) and the other including the biotic data (DNR). Some data were transformed to meet the assumption of normality. We tested whether the variables in a data set (abiotic or biotic) yielded significant discrimination between lakes inhabited by D. parvula versus D. retrocurva and the correlation of these variables with the discriminant function was used to determine the relative contribution of each variable to the group separation. In the abiotic DA we retrieved three variables that significantly discriminated between lakes inhabited by D. parvula or D. retrocurva (lake volume, area, and depth). Since all three variables were mutually correlated, we ran an Analysis of Variance (ANOVA) on lake volume since it had the highest correlation with the discriminant function, with species as the fixed effect. For the invertebrate predators, an ANOVA was used to test for differences between individual predator means indicated by the DA that best discriminate between habitats with D. parvula or D. retrocurva.
We tested whether D. parvula and D. retrocurva differ in their vulnerabilities to habitat-specific predators and if they exhibit morphological traits that are known to affect vulnerabilities towards the respective predators. One experiment used an invertebrate predator (Leptodora kindti) while the other a vertebrate predator (emerald shiners, Notropis atherinoides). In addition, morphological measurements were taken of relative body size and transparency for populations used in the predation experiments.
During late summer (August and September) we collected D. parvula from Buffalo, NY (43° 01' 39.37'' N, 79° 48' 72.39'' W) and D. retrocurva from embayments of western Lake Ontario (43° 11' 22.04'' N, 77° 31' 75.09'' W from Irondequoit Bay, NY in August and 43° 17' 59.83'' N, 79° 48' 16.81'' W from Burlington, Ontario in September). Daphnia was reared for two weeks prior to the experiment at constant densities in 50-ml vials filled with synthetic lake water (96 mg NaHCO3, 60 mg CaSO4 2H2O, 60 mg MgSO4 and 4 mg KCL in 1 L of double distilled water; EPA 2002) with 10 Daphnia per vial at 21°C with a 24 h light photoperiod. All vials were treated with cetyl alcohol (CH3 (CH2)15OH) to prevent the Daphnia from getting caught in the surface tension. Every three days, 1000 μl of Selenastrum capricornutum suspensions were added to each vial of Daphnia as a food source.
Leptodora kindti was field-collected from Lake Ontario (43° 11' 22.04'' N, 77° 31' 75.09'' W), one week prior to the experiment. Leptodora were individually maintained in 50 ml vials filled with synthetic freshwater at 21°C with 24 h light photoperiod. Every three days L. kindti was given zooplankton collected from the field as a food source.
Emerald shiners (from 50 to 64 mm long in standard length) were obtained from a bait shop and were previously field-collected from the Niagara River, Buffalo, NY. Although the source of these predators were from the Niagara River, emerald shiners are often the dominant planktivorous fish in lakes with Daphnia representing one of their main food sources [93, 94]. Prior to the experiment, the fish were kept in 10 L tanks with treated tap water at room temperature (21°C) and fed a mix of freshwater fish food (Tetrafin) and zooplankton collected from the field. Both invertebrate and vertebrate predators were starved for 24 hours prior to the experiments.
For the invertebrate predator experiment, all experimental trials took place in 90 ml plastic cups with 50 ml of synthetic fresh water treated with cetyl alcohol 24 hours prior to the experiment. On the day of the experiment, 20 adult female D. parvula (Buffalo population collected in August) or D. retrocurva (Irondequoit Bay population) were added to a cup, which was then randomly assigned 1 or 0 adult L. kindti (predator and control, respectively). There were 20 replicates of each prey treatment within the predator treatment and 4 replicates of each prey treatment within the control treatment. After 24 hours (21°C, 24 h light photoperiod) the predator was removed and surviving prey were counted.
For the vertebrate predator experiment, all experimental trials took place in 3.8 L plastic containers with 3 L of synthetic lake water at room temperature (21°C). On the day of the experiment, 50 adult females of either D. parvula (Buffalo population collected in September) or D. retrocurva (Burlington population) were added to a container that was randomly assigned 1 or 0 fish (predator and control, respectively). There were 3-4 replicates of the predator treatments and 2 replicates of the control treatments for both prey species on each of 3 days (total: predator 10 replicates, 6 control replicates, for each prey species). The number of replicates in the vertebrate predator experiment was limited by the number of available Daphnia and each replicate used a novel predator. After one hour the predators were removed and the surviving prey were counted.
Prey mortality was equated with the proportion of prey missing. The proportion of prey missing was arcsine square root transformed and analyzed by factorial ANOVA with predator, prey, and the interaction as fixed effects. Post hoc tests with Tukey correction were performed to detect any pairwise differences among mean proportion prey missing between species.
where P I is the initial concentration of prey per experimental unit, P T is the concentration of prey per experimental unit at the end of the trial, X is the number of predators per experimental unit, and T is the duration (in hours) of the trial . Although this equation does not take into consideration reproduction or intrinsic mortality of both prey and predation, in our experiments no reproduction was observed and prey mortality was quantified by the number of prey missing rather than number of dead prey. Although we observed a small proportion of prey missing in the controls treatments, this random effect was likely constant among both control and predator treatment as the higher proportion of prey missing was consistent in the predator treatment replicates.
Following the experiments, morphological measurements were taken on 25 adult females of D. parvula and D. retrocurva for each predation experiment (50 total for each species) that were randomly sampled from the laboratory source populations that provided the Daphnia for each predation experiment (individuals from the populations not used in the experiments). From populations from both predation experiments, we measured the body length (from the top of the compound eye to the base of the carapace), and the total body length (from the highest point of the helmet or head to the tip of the tail spine). Measurements were taken with a dissecting microscope and ImageJ (ver. 1.37, 2006). Since D. parvula cannot produce a helmet, these two measurements were used to quantify the body size of D. parvula relative to that of D. retrocurva, with and without the cyclomorphic features exhibited (helmets and tail spines). We ran separate nested ANOVA's to test for differences in body length and total body length with species as the fixed effect (D. parvula or D. retrocurva), and source population nested within species (locality and time collected for either for the invertebrate or vertebrate predation experiment) as a random effect. Both variables were log transformed to meet the assumptions of normality and homogeneity of variances, although raw means are illustrated in the figures.
From populations of D. parvula and D. retrocurva used for the vertebrate predation experiment, we additionally quantified the relative transparency using a dissecting microscope and the "Histogram" analysis of ImageJ. Each individual Daphnia was photographed in grey scale under identical magnification and conditions, using only transmitted light. The relative transparency was quantified from the mean pixel intensity (0 to 255 for pixel shades from black to white) for the specimen image. Following confirmation that the photo conditions were standardized across both species, an ANOVA was run to test for differences in transparency between D. parvula and D. retrocurva. An identical analysis was also run on transparency only among non-gravid females from both species, as the presence of parthenogenic eggs can affect detection of prey by vertebrate predators . All dependent variables were log transformed to meet the assumptions of normality and homogeneity of variances, yet raw means are illustrated in the figures.
Specimens of D. parvula and D. retrocurva were collected from populations (n = 29) throughout the continental United States and Canada (see Additional file 1, Table S1). Most populations were from areas in Northern United States and southern Canada where the two species' distributions overlap (Fig. 2, Additional file 1: Table S1). We obtained from 1-10 individuals per population (average 5.44 individuals/population).
Total genomic DNA was extracted using Quick Extract (Epicentre). Samples were homogenized in 35-45 μl of Quick Extract solution, then incubated at 65°C for 2 h and 98°C for 10 minutes, and stored at -20°C. One mitochondrial and three nuclear markers were sequenced from the extracted DNA. A 631 bp fragment of mitochondrial protein-coding NADH-2 (ND2) was amplified with the primers (5' - GTTCATGCCCCATTTATAGGTTA - 3') and (5'- GAAGGTTTTTAGTTTAGTTAACTTAAAATTCT-3'). The primers (5'-TTACGAGTCCAGATGGGCTT-3') and (5'- ATCCGTTATGAATCCCTGACTGA-3') were used to amplify a 669 bp fragment of protein-coding HSP90 gene. A 433 bp fragment of the nuclear rab GTPase (F6F12) gene was amplified with the primers (5'- CGTTTCGAATTGGCTTACTGA-3') and (5'- CATCGTTATCTGTCTACGTCTTGAA-3'), while a 534 bp fragment of the translation initiation factor (G6G12) gene was amplified with the primers (5'- AGAAATTCAACATGCCCAAGA-3') and (5'- CGTCGACGAAGTTGACAGTAT-3') . Each Polymerase Chain Reaction (PCR) was 50-54 μl in total and consisted of 4-8 μl of extracted DNA, 5 μl of10 × PCR buffer [50 mM KCl, 1.5 mg MgCl2, 10 mM Tris-HCl pH 8.3, 0.01% (w/v) gelatin], 1.5 μl of DNTP's (2 mM of each), 1 μl of 10 μM of each primer and 1 μl of standard Taq DNA polymerase. Each PCR was conducted on a MJ Thermocycler with the following conditions for the mitochondrial (ND2) gene: 40 cycles at 94°C for 30 s, 48°C for 30 s, and 72°C for 1 min, and a final extension at 72°C for 7 min. The PCR temperature profiles for the nuclear (HSP90) gene were 40 cycles at 94°C for 30 s, 50°C for 30 s and 72°C for 1 min. with a final extension at 72°C for 5 min. For the nuclear (F6F12) gene the PCR temperature profile was 40 cycles at 94°C for 30 s, 58°C for 30 s and 72°C for 1 min. with a final extension at 72°C for 5 min, identical temperature profiles were used for the nuclear (G6G12) gene with the exception of the annealing temperature set at 53°C rather than 58°C. Sequences of all nuclear and mitochondrial PCR products were obtained in both directions by Genaissance Pharmaceuticals or the University of Washington and were both assembled and edited with SEQUENCHER ver. 4.2 (Gene Code) then manually aligned in SE-AL 2.0 .
An individual was considered heterozygous when overlapping peaks (lower peak > 95% of higher peak) were observed at any given site on the sequence electropherogram. Individuals with multiple heterozygous sites were cloned with the Invitrogen TOPO TA kit to determine the different alleles at each site with observed heterozygosity. The QIAprep Spin Miniprep Kit (QIAGEN) was used for plasmid purification and the respective primers for each particular locus were used to sequence the cloned inserts. For each individual, six cloned fragments were sequenced to help detect cloning artifacts . The combined data for all four makers was a 2268 bp sequence with both alleles represented for each individual for the nuclear genes. All individuals sequenced for the F6F12 and G6G12 markers were also represented with sequences for the HSP90 and ND2 markers. DNA quality appeared to prohibit sequencing of all genes for all individuals. Therefore, we conducted analyses for 1) individuals represented by sequences for only the mitochondrial ND2 gene, 2) individuals represented by sequences with only the three nuclear genes, and 3) individuals represented by sequences for all four genes (total evidence).
The number of individuals sequenced from each marker used in the analyses can be seen in Additional file 1: Table S1. Each individual was represented by two copies of each sequence to account for both alleles at a heterozygote site in the nuclear markers, while the mitochondrial (ND2) gene was replicated twice to correspond with the two copies of each nuclear gene sequence for the analyses. We performed a phylogenetic analysis on the total evidence data set (all four genes combined) using Maximum Likelihood (ML) method with RAxML  and a file to partition the alignment by codon position and gene. For all analyses a GTR substitution model was implemented for the analysis and unique model parameters were optimized and applied for each partition. Support values were estimated using 100 bootstrap replicates.
Because the transition from polyphyly to monophyly is continuous, a categorical description of divergence can fail to identify significant divergence before monophyly is reached . Thus, detecting divergence in closely related or recently diverged species phylogenetically can be difficult if monophyly has not been reached. The recently described Genealogical Sorting Index (gsi) allows one to detect significant genealogical divergence before monophyly by quantifying the relative degree of exclusive ancestry of a labelled group on a rooted tree topology. Since branch tips along a tree generally represent gene copies of a given locus, the value of gsi is the degree of genealogical exclusivity among the sampled gene copies. The calculated gsi value ranges from 0-1, with 0 representing panmixia and 1 representing monophyly and permutation tests provide a statistical test of significance for this value .
We calculated the gsi value for each species (D. parvula and D. retrocurva) and tested the null hypothesis that the gene copies labelled D. retrocurva and D. parvula form a single intermixed group. Because strongly unbalanced group representation can result in a decreased power to detect significant exclusive relationships  we determined the gsi for the mitochondrial data, the nuclear data, and the total evidence data after randomly removing D. retrocurva individuals from the data to achieve a balanced sampling scheme, since this species was over represented compared to D. parvula.
A median-joining haplotype network was constructed for the mitochondrial data and the three combined nuclear gene sequence data using NETWORK ver. 18.104.22.168 . In order to further estimate divergence time we used the isolation with migration- analytic (IMa) model, which was designed to analyze recently separated populations (or in this case, species) not at equilibrium. This model estimates six demographic parameters including gene flow rates and time of divergence while generating relative likelihoods/posterior probabilities . Unfortunately we could not get reliable migration or gene flow estimates in IMa likely due to the increased number of parameters in the model (J. Hey communication); therefore, we used IMa solely for estimating divergence time. Since the degree of current gene flow between D. parvula and D. retrocurva is unknown we ran two IMa analyses, one with the migration prior set to zero migrants per generation and another with ten migrants per generation between the two species which we treated as populations. For both analyses, we used the Hasegawa-Kishino-Yano model  and the Infinite Sites mutation model for the mitochondrial and nuclear markers, respectively. The mutation rates entered in the model were the number of bp in the particular marker multiplied by 3.3 × 10-6 mutations/site/year for the mitochondrial marker and 1.05 × 10-7 mutations/site/year for the nuclear markers , assuming 5 generations per year; previous estimates adjusted to seasonal duration of localities [106, 107]. The analysis included 30 million generations post burn-in (1 million generations) with six Metropolis-coupled Markov chains with the first heating parameter for the linear heating scheme set to 0.05. Tajima's D test and the four-gamete test (DNASP ver. 4.0) indicate all markers are consistent with the expectations of neutrality and the nuclear markers show no evidence of recombination, both assumptions of the IMa model.
Sequences from the rapidly evolving mitochondrial ND2 gene were used to determine the demographic history of both D. parvula and D. retrocurva lineages. Three mismatch distributions were calculated from the number of differences between pairs of haplotypes within and between species in ARLEQUIN . The Sum of Square Deviations (SSD) was computed for each species separately then tested whether the observed distributions deviated significantly from those expected under the population expansion model using 10000 permutation replicates. We used the parameter τ to estimate the time since the expansion (t) using the equation t = τ/2 u, where u is the mutation rate per sequence per generation  assuming a mutation rate of 6.6 × 10-8/site/generation  and five generations/year.
We ran hierarchical Analysis of Molecular Variance (AMOVA) in ARLEQUIN  to determine the total genetic variation explained by differences (i) between the two species D. parvula and D. retrocurva, (FCT), (ii) among populations within a species, (FSC), and (iii) within a population (FIS).
We thank Howard Lasker, Charles Mitchell, Howard Reissen, and Steven Juliano for their comments and suggestions with analyses and interpretation. We also thank Dr. Kenton Stewart, Laura Hovind, Seiji Ishida, and Angela Omilian for assistance with specimen collections, data collection and analyses. This research was supported by an NSF grant (DEB- 0331095) to D. J. T..
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.