Ethics statement
Local administration has been contacted at the beginning of the study. Collection of zooplankton (Daphnia) did not require specific permissions as these samples were obtained from unprotected lakes that are open for public activities. Our study did not involve the use or collection of endangered or protected species.
Study sites and field sampling
Daphnia individuals were sampled from ten flooded gravel pits in and around Munich (Germany; Additional file 1: Table S1). The lakes were created between 1930 and 2000, most are small and shallow, and the maximum distance between them is about 30 km. Geographical coordinates, lake age, lake size, maximum depth, total phosphorus concentration (Ptot), nitrogen concentration (N) and Daphnia densities (averaged across monthly densities from April to November 2011) are given in the supplementary material (Additional file 2: Table S2). The phosphorus and nitrogen data were provided by Andreas Scholz from the Wasserwirtschaftsamt München (http://www.wwa-m.bayern.de/). These measurements were taken in spring, from the surface to the bottom of the lake at one-meter intervals (average Ptot and N values were then calculated). Zooplankton samples were collected monthly from April to November 2011 (eight samples per lake; sampling interval was about 30 days). A 95-μm plankton net was hauled through the whole water column at two different sites per lake within the deep basin (depth was measured using a portable depth sounder); the two samples were then pooled and preserved in 95% ethanol. A random subsample was used for density counts (only adult females were taken into consideration). Then, ca. 50 adult females from the D. longispina complex were chosen randomly per sample for genotyping. Since Daphnia densities were sometimes very low, a complete time series of eight samples could be genotyped for only three lakes (FASA, LERC and WALD; for lake abbreviations see Additional file 1: Table S1). Two to seven samples were genotyped for the seven remaining lakes, resulting in 1934 Daphnia individuals being genotyped in total.
Microsatellite genotyping
The DNA of each individual Daphnia was extracted and genotyped at 15 microsatellite markers [37], in two multiplex polymerase chain reactions, following a protocol described elsewhere [18]. PCR products were analysed on an ABI PRISM 3700 capillary sequencer using a LIZ 500 labelled size standard. Genotypes were scored using GeneMapper version 3.7 (Applied Biosystems). Alleles at each locus were defined by their fragment length (in base pairs). Across different runs of genotyping, the consistency of alleles was checked, with locus-specific patterns of one reference genotype used in each run. This enabled us to adjust alleles with small differences in fragment length among different runs. In addition, there was no evidence of scoring errors due to stuttering, large allele dropout, or the presence of null alleles [38], as indicated by respective tests using MICRO-CHECKER 2.2.3 (104 permutations, [39]).
Taxon assignment
To display the genetic relatedness of multilocus genotypes (MLGs), factorial correspondence analysis (FCA) was applied in GENETIX 4.05 [40], on all unique MLGs together with 49 well-defined reference genotypes (same set of genotypes as in [18]), representing three parental species (D. cucullata, D. galeata, D. longispina) and their interspecific hybrids. It was confirmed that only D. galeata and D. longispina parental species were present in our dataset. Bayesian statistics were then used (in NewHybrids 1.1, [20]) to assign individuals to one of six possible classes (i.e., two parental species, F1 and F2 hybrids and both backcrosses). The probability threshold for assignments was set to 95% (106 iterations after a burn-in of length 106). All the following calculations of taxon-specific parameters are based on the taxon assignment in NewHybrids.
Spatial patterns in taxon composition
To test if the taxonomic similarity of Daphnia communities corresponds with the geographical proximity of the lakes, the correlation between pairwise Euclidean distances based on the taxon composition and respective pairwise geographical distances between sampling sites was evaluated using a Mantel test (5000 permutations, Past v2.07, [41]). Two estimates of Euclidean distances were derived in Past v2.07 for the spring samples (i.e., April or May), based on either the frequency, or the presence/absence, of five Daphnia groups. These groups included all classes as identified by NewHybrids (i.e., D. galeata, D. longispina, F1 hybrids and F2 hybrids; backcrosses were not detected) and an additional group of individuals that could not be identified by this method. All ten lakes were included in these tests.
Environmental preferences of parental taxa
A principal components analysis (PCA) of five environmental variables (i.e., lake age, lake size, maximum depth, total phosphorus concentration (Ptot), nitrogen concentration (N)) and mean Daphnia density (see Additional file 2: Table S2) was used to relate Daphnia species occurrence to combinations of these parameters. Variables that did not conform to normality were log-transformed. To test if Daphnia species were nonrandomly distributed along each component of the PCA, a one-way ANOVA was used with the presence/absence of each Daphnia species (i.e., D. galeata or D. longispina) as a main factor. The species presence/absence category was assigned based on the whole-year dataset (Additional file 1: Table S1). All ten lakes were included in these tests.
Temporal changes in taxon composition
The relative frequencies of five Daphnia groups (i.e., four classes: D. galeata, D. longispina, F1 hybrids and F2 hybrids identified by NewHybrids, and one group of individuals not identified by this method) were compared between spring (April or May) and autumn (October or November or, in the case of FELD, September), by applying a Monte Carlo approach with 105 simulation runs [42]. Since separate tests were run for each of the seven lakes tested here (BOHM, FASA, FELD, HEIM, LERC, LUSS and WALD), sequential Bonferroni corrections [43] were applied when interpreting the results.
Genetic distances among spatially and temporally isolated populations of parental species and hybrids
Pairwise Nei’s genetic distances [21] were calculated in GENALEX 6 [44], among populations belonging to D. galeata, D. longispina or hybrids (F1 and F2 NewHybrids classes were pooled into one group of “hybrids”, because of the low sample size of the F2 class). Population samples originating from different months were considered separately, and only populations with sample size N ≥ 6 were included. Then, genetic similarities among population samples were displayed using the Unweighted Pair-Group Method with Arithmetic Mean (UPGMA), as calculated in MEGA 4 [45].
Gene flow and local production of hybrids
The number of migrants between lakes was calculated per generation (Nm) to estimate gene flow among populations of the two parental species (i.e., D. galeata and D. longispina) and the hybrids (F1 and F2 NewHybrids classes were pooled into one group of “hybrids”, because of the low sample size of the F2 class), respectively (in GENETIX 4.05). The whole-year dataset was pooled per taxon and lake. For each of the three taxa, only lake samples with N > 10 were considered to represent a local population and included in this analysis. Then, one-way ANOVA and Post-Hoc Tests were used to test for taxon-specific differences in the number of migrants (in SPSS 20.0).
In addition, it has been tested if hybrids are more likely formed by the parental species inhabiting their home lake (i.e., local production of hybrids) than by the parental species from other lakes. Thus, for each individual hybrid genotype, the log-likelihood with which it was assigned to parental groups (specific parental species from the home lake versus pooled populations of the same parental species from the other lakes) was calculated in Arlequin 3.0 [46]. To distinguish statistical differences in respective log-likelihoods of hybrid assignments, a paired t-test (103 bootstrap replications) was run in SPSS 20.0. Although hybrids and respective parental species co-existed in FELD and LUSS, the test above could only be performed for LUSS hybrids, because only one hybrid genotype was detected in FELD.
Test for clonality
To test the null hypothesis that an MLG encountered more than once was the result of sexual recombination, instead of clonal propagation, the Psex value was calculated in GENCLONE 2.0 [47]. The calculations were performed individually for each abundant taxon and lake. Only individuals characterized at all 15 microsatellite loci (or with missing data at a single locus, SwiD2; many D. longispina individuals could not be amplified at this locus) were included in these and the follow-up analyses of clonal structure (i.e., 1856 of 1934 genotyped Daphnia, see Additional file 1: Table S1).
Temporal changes in clonal composition
To test if differences in clonal composition increase with the time interval between samples, the clonal turnover rate was calculated between sampling dates (for all possible pairwise comparisons), per taxon and lake, by computing complements of the Morisita-Horn index (1-MH [48]; see also [49]) using the program SPADE [50] and 104 bootstrap replications. This turnover rate takes into account not only changes in the presence of individual MLGs but also changes in MLG frequencies. Turnover rate varies from 0 (complete similarity) to 1 (no similarity). Then, the correlation between the pairwise turnover rate in clonal composition and the corresponding pairwise time differences between sampling dates (calculated in monthly intervals) was evaluated by a Mantel test (5000 permutations, Past v2.07). Finally, evidence was sought for the existence of an overall decrease in clonal diversity during the course of the growing season (i.e., clonal erosion). Clonal diversity was measured per taxon and sampling date as: (i) clonal richness R, and (ii) the inverse of the Simpson index. Clonal richness R was calculated as R = (G-1) / (N-1), where G is the number of genotypes and N represents sample size [51]. R varies from 0 (all individuals belong to one clone) to 1 (all individuals are different). The inverse of the Simpson index was calculated as 1/(∑Pi2), where Pi is the proportion of the ith clone. This index thus includes additional information about the distribution of clones within samples. The minimum value of this index is 1 (all individuals belong to one clone) whereas the maximum value is the number of unique clones in the sample. Linear regression analyses were applied to test for the effect of time on clonal diversity (i.e., either clonal richness R or the inverse of the Simpson index) using R-software [52]. The tests were run only on those populations for which samples spanning the entire growing season were available: D. galeata from FASA, D. galeata from LERC and D. longispina from WALD (a minimum sample size was set of 30 individuals).