Aquatic insects are widely used as indicator taxa to assess the ecological state of streams and to evaluate the success of stream restoration projects. Information on intraspecific genetic diversity and population connectivity is often lacking for such indicator taxa. However, these parameters are of critical importance for restoration plans and conservation management because: 1) species sometimes consist of several cryptic species and 2) species can recolonize only those restored habitats within a reachable distance from their source populations. Gene flow generally cannot be observed directly, and molecular markers provide a reasonable alternative to assess the dispersal potential and evaluate species' genetic diversity. We investigated the genetic diversity and dispersal potential of the predatory stonefly Dinocras cephalotes using 323 specimens from 29 populations in the Sauerland, a low mountain range in Germany. We used a 658 base pair (bp) fragment of the mitochondrial cytochrome c oxidase subunit I gene (COI) and found 2 distinct and diverse haplotype groups, which were shared across most populations. The groups were separated by a minimum intraspecific P-distance of 4.3%, suggesting historic isolation and possible presence of cryptic species. However, complementing analyses of the nuclear Wingless gene and 3 newly developed microsatellite markers clearly showed that individuals from both COI haplotype groups are interbreeding, and therefore, D. cephalotes is considered a single valid species. Population comparisons indicated high connectivity among all populations, with only a few individual populations showing signatures of isolation. Based on the molecular data, we conclude that dispersal is primarily achieved by the adult females of D. cephalotes.
Human activity has dramatically altered and degraded stream ecosystems (Poff et al. 2007) and poses a threat to global freshwater biodiversity (Vörösmarty et al. 2010). The loss of biodiversity threatens ecosystem functioning (Vaughn 2010), with potential direct negative consequences for the provision of ecosystem services (Cardinale 2011). The Water Framework Directive (WFD) of the European Union explicitly obliges its member states to counteract degradation of freshwater ecosystems and demands good ecological and chemical status of surface waters by 2015 (Directive 2000/60/EC, Annex V; European Union 2000). In a recent report on the status of Europe's waters, the European Environment Agency concluded that most surface-water bodies in Europe are unlikely to meet this target by 2015 (EEA 2012), and many management and restoration activities have been launched to implement this ambitious Directive.
The criteria for attaining good ecological status have not yet been fully developed, but the recovery of intact freshwater communities is the key to improve ecological conditions in lotic systems (Palmer et al. 1997, 2010, Jähnig et al. 2009, Feld et al. 2011). A primary assessment tool for quantifying the ecological status of streams is analysis of the biodiversity and abundance of biological indicator species, in particular macroinvertebrates (Hellawell 1986, Metcalfe 1989, Rosenberg and Resh 1993, Hering et al. 2006). Once restoration projects are completed, the native bioindicator organisms must be capable of reaching the restored ecosystems from source populations (Palmer et al. 1997, Lake et al. 2007, Smith et al. 2009). Successful recolonization depends on several additional variables, such as the individual species' life cycle (holo- vs merolimnic) and its duration, physical dispersal traits (flying, crawling, active, and passive), behavioral patterns, the spatial distribution of source populations in the stream network, and the presence of barriers that hinder dispersal (Smith et al. 2009). Macroinvertebrate species often may be unable to re-colonize restored habitats successfully, despite favorable habitat conditions at the restored site (Lake et al. 2007). In such cases, the absence of macroinvertebrate indicator taxa is not automatically indicative of poor habitat quality and misleads stream-quality assessment.
Empirical or experimental data on the dispersal abilities of aquatic insects are scarce and, therefore, often neglected in restoration plans (Palmer et al. 1997, Smith et al. 2009). Studies of aquatic insects with stable isotopes (Coutant 1982, Briers et al. 2004), mark-and-recapture experiments (Stettmer 1996, Has-sall and Thompson 2011), and light traps (e.g., Kovats et al. 1996) indicate that long-distance dispersal or passive drift >1 km from the stream channel is possible for winged insects. Although a strong preference of insects to stay in the vicinity of the stream channel can be observed (90% of adult stonefiies stayed within the stream channel; Briers et al. 2002), rare but successful long-distance flights may connect populations. In the absence of direct observations, molecular tools allow assessment of the connectivity of populations with comparisons of allele diversity and frequencies among populations (Hughes et al. 2008). These markers also help identify overlooked or cryptic species, which show no or only subtle morphological differences (Pfrender et al. 2010, Zhou et al. 2010).
We analyzed partitioning of genetic diversity in the predatory stonefly Dinocras cephalotes (Curtis 1827) in the Sauerland region, a low mountain range in western Germany. The larvae of D. cephalotes inhabit cold and fast-flowing streams and are reliable indicators of good water quality (Eiseler and Enting 2012). Dinocras cephalotes has a life cycle of ∼3 y (Iannilli et al. 2002), and only the female imagines have fully developed wings (Tierno de Figueroa et al. 2006). The dispersal potential of D. cephalotes at local and regional scales has been questioned because they were described as clumsy flyers (Ketmaier et al. 2001) and population subdivisions were reported even on very small geographical scales (Ketmaier et al. 2001).
The Sauerland region has been anthropogenically influenced since medieval times by agriculture, forestry, ore mining, and metal production, which in concert, led to severely impacted stream networks. However, many headwater streams in the Sauerland region are chemically and ecologically classified as being in good condition. Isolated and genetically depauperate headwater populations of D. cephalotes might be expected because of the heavy influence of several anthropogenic stressors (e.g., pollution and fragmentation by weirs, dams, and hydromorphological alterations) on higher-order streams in the network. This isolation may be especially strong for D. cephalotes because flight capability is restricted to females.
The aim of our study was to use mitochondrial and nuclear deoxyribonucleic acid (DNA) markers to test whether populations of D. cephalotes in different headwater streams are isolated and genetically depauperate, which in turn, would increase demographic stochasticity and the risk of local extinctions. The regional dispersal potential, which is the prerequisite to recolonize restored stream habitats, of D. cephalotes and its value as an indicator taxon for assessing restoration measures were assessed.
Dinocras cephalotes populations in headwater streams were sampled, mainly in the Ruhr river basin of North Rhine-Westphalia, Germany (rivers Ruhr and Lenne; Fig. 1, Table S1 (TableS1.pdf)) from May to June 2010–2012 and stored in 80% ethanol at -20°C. A total of 323 specimens from 29 populations were analyzed with molecular methods.
Microsatellites were developed from an unen-riched sequence library that was sequenced on 454 GS Junior sequencer (Roche, Basel, Switzerland). One microgram of high-quality DNA was used for library creation according to the manufacturer's protocol. Resulting reads were quality-controlled (FastQC, version 0.10.1; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), analyzed for potential contamination (Blast+, version 2.2.26; Camacho et al. 2009), and subsequently assembled using MIRA (version 188.8.131.52; Chevreux et al. 1999) to exclude potential multicopy sequences (settings as in Leese et al. 2012 except -AL: mrs=70:mo=20 -CL:pec=off ). Microsatellites were identified with Phobos (version 1.0.6; Mayer 2006) and a custom R script (version 2.15.1; R Project for Statistical Computing, Vienna, Austria) was used to select the best candidate sequences based on flanking region length, repeat motif, and number. Sequences were inspected in Geneious Pro (version 6.0.5; Kearse et al. 2012). Primers were developed with the Primer3 plug-in (Rozen and Skaletsky 2000).
All specimens were identified and photographed prior to DNA extraction. DNA was extracted following a modified salt-extraction protocol (Sunnucks and Hales 1996). Negative controls were included for DNA extractions and polymerase chain reactions (PCR). A fragment of the mitochondrial barcoding gene cytochrome c oxidase subunit I (COI) was amplified with standard invertebrate primers HC02198 and LCO1490 (Folmer et al. 1994) in a reaction consisting of l × PCR buffer, 2.5 mM MgCl2, 0.2 mM deoxyribonucleo-tide triphosphates (dNTPs), 0.5 (µM of each primer, 0.02 U/(µL Euro Taq (EuroClone, Milano, Italy), 1 (µL DNA, filled up to a total volume of 25 (µL with high-performance liquid chromatography (HPLC) H2O (PCR program: 94°C/120 s, 36 cycles of [94°C/20 s, 46°C/30 s, 72°C/60 s], 72°C/7 min). PCR success was confirmed with agarose gel electrophoresis, and samples that could not be amplified were repeated using HotMaster Taq (5Prime; Gaithersburg, Maryland; parameters identical as for Euro Taq, but with 65°C extension temperature). A fragment of the nuclear Wingless gene was amplified from 68 samples of 10 populations (illustrated with an asterisk in the haplotype map; Fig. 1) according to the protocol described by Pauls et al. (2008). Populations for additional Wingless analysis were selected to cover a wide geographic range and, if possible, to include individuals belonging to both COI haplotype groups. Ten µL of PCR product were purified enzymatically with 0.5 (µL Exonuclease I (20 U/(µL; Thermo Fisher Scientific, Waltham, Massachusetts) and 1 (µL FastAP (1 U/µL, Thermo Fisher Scientific) by incubating in a thermocycler at 37°C for 15 min followed by 96°C for 15 min prior to sequencing (Werle et al. 1994). Bidirectional Sanger sequencing was carried out on an ABI 3730 sequencer (Applied Biosystems, Carlsbad, California) by GATC Biotech (Kon-stanz, Germany). Ambiguous and short sequences were resequenced.
Microsatellite markers were optimized with a gradient PCR (annealing temperature range: 45–69°C) and PCR enhancers (dimethylsulfoxide [DMSO]: Carl Roth, Karlsruhe, Germany; Betaine: Sigma-Aldrich, Steinheim, Germany). The PCR settings were as described above for COI, but with a concentration of 0.2 µM for each forward and reverse primer, 0.05 [µM of the tailed M13 primer (5′-CACGACGTTGTAAAA CGAC-3′), and 0.02 U/µL Euro Taq (PCR program: 94°C/120 s, 36 cycles of [94°C/20 s, 51–63°C/30 s, 72°C/60 s], 72°C/45 min). Primers that amplified reliably for a subset of samples were used for all samples. Allele sizes were determined by acrylamide gel electrophoresis on a Li-Cor analyzer 4300 with the software Saga2 GT (Li-Cor Biosciences, Lincoln, Nebraska). Alleles that could not be identified reliably were rerun or scored as missing data.
Sequence data analysis
Sequences with sufficient length and quality were assembled and an alignment was constructed using the MAFFT plugin (version 1.3; Katoh et al. 2002) for Geneious. Uncorrected genetic p-distances between haplotypes were calculated with MEGA (version 5.05; Tamura et al. 2011). To test for genetic fixation between populations, pairwise FST and ΦST estimators were calculated for the COI data set using Arlequin (version 184.108.40.206; Excoffier et al. 2005).
Statistical significance was assessed with 1000 iterations and significance level was Bonferroni corrected. Isolation by distance was tested with a Mantel test (10,000 replications; R package ade4; Dray and Dufour 2007) using the genetic differentiation measures FST and ΦST with 3 distance measures: direct distance between populations, shortest distance following the streams, and elevation differences between populations (calculated with QGIS, version 1.8; Quantum GIS Development Team; http://qgis.osgeo.org). Population partitioning with east/west and grouping by catchments was tested with an analysis of molecular variance (AMOVA) as implemented in Arlequin. Grouping was used to test for local adaptations resulting from effects related to altitude and isolation between catchments in case of poor dispersal. A minimum spanning network was calculated with Arlequin and visualized with HapStar (version 0.6; Teacher and Griffiths 2010).
Wingless sequences were assembled like the COI sequences, and a minimum spanning network generated as described above. Wingless haplotypes were compared to the respective COI haplotypes of the 63 tested specimens. The Wingless marker is a nuclear gene and was sequenced to validate the patterns found with the mitochondrial marker COI. A relatively small sample size for the Wingless gene was sufficient to test whether the patterns of both markers were similar.
Microsatellite data analysis
The microsatellite data were checked for scoring errors with MicroChecker (version 2.2.3; Van Ooster-hout et al. 2004) and for deviations from Hardy-Weinberg and linkage-equilibrium with Arlequin. To measure genetic differences between populations FST values were calculated with GenoDive (version 1.0b23; Meirmans and Van Tienderen 2004). Destwas calculated with the R package DEMEtics (version 0.8–5; Jueterbock et al. 2012) with Bonferroni-corrected p-values (populations with a sample size = 1 were excluded from the analysis). FST between the 2 main haplotype groups for the COI marker was calculated by creating 2 artificial populations, each containing the microsatellite data for all individuals of haplotype group A or group B (GenoDive, 1000 iterations). The Mantel test and AMOVA were calculated as described for the COI data. In addition, the microsatellite data were analyzed for population clustering using STRUCTURE (version 3.2.4; Pritchard et al. 2000, Falush et al. 2003; default settings, burn-in = 10,000 followed by 50,000 Markov Chain Monte Carlo steps) and the most likely number of clusters was determined with the Evanno method (Evanno et al. 2005) as implemented in STRUCTURE HARVESTER (version 0.6.93; Earl and vonHoldt 2011; http://taylor0.biology.ucla.edu/structureHarvester/ ). Last, a principal component analysis (PCA) was carried out in the R package Adegenet (Jombart 2008) for data conversion (function scaleGen and dudi.pca) as implemented in the package ade4.
Haplotype groups (COI and Wingless)
Reliable COI sequences of 658 bp length were obtained for 307 specimens (GenBank accession numbers KF410897-KF410943). The other sequences were excluded because of poor read quality, short read length, or double peaks. Some sequences that still showed double peaks after PCR and sequencing reactions were repeated and may hint at the presence of pseudogenes (numts) or heteroplasmy. Reliable sequences of 400 bp length were obtained for all of the 68 samples analyzed for the Wingless gene (GenBank accession numbers KF442621-KF442625). However, 2 single-nucleotide polymorphisms (SNPs) were observed in 5 sequences. The individual peaks in the chromatograms showed very exact overlap, so phasing was not possible, and the 5 sequences were discarded from the data set.
A total of 47 unique COI haplotypes were found. These clustered into 2 diverse groups separated by p-distances of 4.3 to 5.2% (Fig. 2A). The Wingless gene showed 1 major haplotype group that was less diverse (only 5 alleles) than the COI marker (Fig. 2B). Most Wingless alleles were shared by individuals from both COI haplotype groups (Fig. 2B). The COI haplotypes had a relatively homogeneous distribution across the study area (Fig. 1). A very weak but significant differentiation existed between the east/west groups (AMOVA, ΦCT = 0.065, p < 0.001; Table 1). Only a few pairwise population comparisons had significant FST values (uncorrected mean FST = 0.1858, σ= 0.0751, p = 0.05, n = 34 of 406; Fig. S1 (FigureS1.pdf)) or significant ΦST values (uncorrected mean ΦST = 0.2960, σ = 0.1605, p = 0.05, n = 32 of 406; Fig. S2 (FigureS2.pdf)). Only 1 pair (FST values for GR and E03) remained significant after Bonferroni correction. The populations E13, E19, GR, NH showed slightly higher differentiation estimates than other populations ( Fig. S1 (FigureS1.pdf), S2 (FigureS2.pdf)). The differentiation values were unreliable for the 3 populations represented by single specimens (RU2, NH, BB) and were not considered further. A weak but significant positive correlation was found between ΦST and differences in altitude between populations (Mantel test, r = 0.1692, p = 0.0097; Fig. S3 (FigureS3.pdf)).
Microsatellite data analysis
Four of 15 developed microsatellite markers were used for all samples (Table 2). Possible stuttering problems were identified for markers C1 and L1 and possible null alleles for the markers C1, C2, and L1 (MicroChecker), but bands of all markers were clearly identifiable despite slight stuttering. All markers deviated slightly from Hardy-Weinberg equilibrium, but only markers C1 and L1 showed very strong deviations (Table 2). L1 was excluded from the data set because of strong heterozygote deficiency and possible linkage to markers C1 and C2.
Results of 2 analyses of molecular variance (AMOVA) for the cytochrome c oxidase subunit I (COI) marker with grouping by catchments and grouping according to geographical position (east/west).
With microsatellite markers, only a few pairwise population comparisons showed significant FST (uncorrected mean FST= 0.009, σ = 0.039, p = 0.05, n = 36 of 325; Fig. S4 (FigureS4.pdf)) or significant Dest values (uncorrected mean Dest = 0.042, σ = 0.091, p = 0.05, n = 50 of 325; Fig. S5 (FigureS5.pdf)). The populations E05, E20 and SO showed slightly higher differentiation values than other populations ( Figs. S4 (FigureS4.pdf), S5 (FigureS5.pdf)) and did not match the populations differentiated in COL Six comparisons had significant Dest values after Bonferroni correction. Most variance was found within populations (AMOVA; Table 3). River distance and Destvalues were weakly but significantly correlated (Mantel test, r = 0.1563, p = 0.0259; Fig. S6 (FigureS6.pdf)). The structure analysis for all populations did not indicate distinct clusters (Fig. 3), and the data analysis with the Evanno method confirmed that one cluster is most likely. The PCA analysis with the microsatellite data clustered all populations together (overlapping ellipsoids; Fig. 4A). Some populations showed deviations from the cluster. However these populations did not show a strong differentiation when analyzed for FST, ΦST, or Dest (COI and microsatellite data). Results of the same PCA analysis also were shown to group by COI haplo-types A and B instead of populations, and both hap-lotype groups showed a clear overlap (Fig. 4B). The microsatellite-based FST between the 2 haplotype groups was not significant (FST = -0.006, p = 0.967) indicating panmixia among members of the 2 groups at the nuclear level.
Overview of developed primer sequences used on 323 Dinocras cephalotes samples. Primers with an M13 extension at the 5′ end are indicated by an asterisk. Hardy-Weinberg equilibrium proportions of expected (He) and observed (H0) heterozygosity were calculated. Significant deviations are indicated by asterisks. * = p < 0.05, *** = p < 0.001. Temp = temperature, F = forward, R = reverse, A = adenine, T = thymine, G = guanine, C = cytosine.
Evidence for cryptic species
An implicit requirement when analyzing the dispersal potential of a target species is that the candidate species does not consist of a complex of cryptic species. An obvious result of the COI analysis of D. cephalotes is the presence of great intraspecific distances (4.3–5.2%) and a prominent barcoding gap between members of 2 haplotype groups (groups A and B). Such genetic signatures often indicate the presence of cryptic or unrecognized species (Hebert et al. 2004). Other plecopteran taxa analyzed so far have shown intraspecific COI distances below and above the typical barcoding gap threshold of 2 to 3% distance (Sweeney et al. 2011, Zhou et al. 2009, 2010), with intraspecific distances of up to 5.8% (Mynott et al. 2011). These studies relied only on COI data, so the possibility that they were dealing with cryptic species could not be ruled out. The variability of D. cephalotes is at the upper margin of values reported as intraspecific distances in the literature (intraspecific uncorrected distance up to 5.2%). A barcode gap of 4.3% between the 2 diverse COI haplotype groups A and B suggests the presence of cryptic species.
Results of 2 analyses of molecular variance (AMOVA) for 3 microsatellite markers with grouping by catchments and grouping according to geographical position (east/west).
However, the barcode gap for the D. cephalotes populations also could be the result of historic isolation and independent lineage sorting in glacial refugia without subsequent reproductive isolation. Therefore, the absence of recombination within mitochondrial genes would lead to the persistence of these historically accumulated differences in secondary contact even under panmixia. Historic isolation has been discussed as a primary force underlying contemporary genetic variation in other aquatic insects (Pauls et al. 2006, Lehrian et al. 2010, Bálint et al. 2011, Alp et al. 2012, Theissinger et al. 2012). The central question to be addressed with regard to such prominent differences in mitochondrial DNA is whether members of these groups still interbreed successfully in secondary contact. If a reproductive barrier had evolved (either because of, e.g., pleiotropic effects in small refugial populations or a specific differential selection regime), the genetic signatures of isolation also should be found in nuclear DNA markers. If interbreeding were still possible, recombination should homogenize the accumulated nuclear differences and lead to differing patterns between mitochondrial and nuclear DNA. To test both hypotheses, the nuclear Wingless marker was sequenced for a subset of individuals from both COI haplotype groups and compared to the COI data. In addition, COI data were compared to allele frequencies of 3 microsatellite markers. Individuals from both haplotype groups shared the same Wingless haplotype and the 3 microsatellites revealed no differences in allele frequencies between members of groups A and B. Thus, the hypothesis that haplotype groups A and B represent genetically distinct cryptic or unrecognized species was rejected.
Connectivity of D. cephalotes populations
The COI haplotypes show a relatively homogeneous distribution across all populations, with a very weak differentiation between the eastern (higher altitude) and western (lower altitude) populations as revealed by significant FCT in the AMOVA. However, only a few pairwise FST and ΦST comparisons showed significant differences even across larger geographical distances. These results indicate connectivity between populations. The small differences between eastern and western populations also could be explained by systematic differences in genotypes related to altitude as shown by the Mantel test.
The results for the 3 microsatellite loci are largely consistent with the results of the COI analysis, indicating stable population connectivity with only a few significant FST and Dest values. No patterns of isolation were revealed by the Bayesian clustering analysis in STRUCTURE or the PCA. However, detailed predictions about the number of migrating individuals per generation cannot be made because of the limited number of microsatellite loci investigated.
A very weak but significant correlation between Dest and river distance was found with the Mantel test. In addition, the differences between eastern and western COI populations found with the marker COI with the AMOVA were not found in the microsatellite data. These results contradict the COI-based results because gene flow among all groups, as indicated by the AMOVA, would affect the distribution of mitochondrial markers. The observed pattern of small but significant mitochondrial differentiation between lower-altitude western and higher-altitude eastern populations might be caused by directional selection acting on the mitochondrial genome (Ballard and Whitlock 2004) or the limited number of micro-satellite loci investigated.
Dispersal strategy of D. cephalotes
Different dispersal strategies are described for aquatic insects (Hughes et al. 2009). Dispersal of D. cephalotes by larval migration within connected streams seems unlikely because a stronger correlation between population differentiation and river distance would then be expected (Hughes et al. 2008). Furthermore, earlier studies based on genetic data showed that most aquatic insects in which both sexes are winged have sufficient dispersal abilities to maintain population connectivity across rivers (Hughes 2007, Hughes et al. 2008).
A study on D. cephalotes with isoenzyme markers found patchy gene flow with some gene flow between the populations of the rivers Aniene, Nera, and Velino in Italy (Ketmaier et al. 2001). Both the Nera and Aniene meet the Tiber River below 100 m elevation, so population connectivity by rivers is unlikely. Therefore, it seems plausible that dispersal is accomplished primarily by flying female imagines. Patchiness in population connectivity might be caused by habitat heterogeneity and land use between river catchments (Smith et al. 2009), as could be the case for D. cephalotes populations in the Sauerland region. Isotope labeling showed that individuals in another stonefly species flew to adjacent catchments at a distance >500 m (Briers et al. 2004).
Dinocras cephalotes has been described as a poor flyer (Ketmaier et al. 2001), but successful migration of only a few ovigerous females per generation may be sufficient to maintain gene flow. COI haplo-types have a relatively homogeneous distribution across the study area and the haplotype groups show a high diversity (47 different haplotypes), so effective population size probably is high because of recurrent exchanges. Thus D. cephalotes cannot be considered as genetically depauperate or endangered in this region. Dinocras cephalotes is a predatory stone-fly (Bo et al. 2007), so its density in a habitat is controlled by prey availability. If populations were indeed small and isolated, they would quickly lose genetic diversity by genetic drift, which, in turn, would be detected quickly with the markers applied in our study.
Two genetically distinct COI haplotype groups were found for D. cephalotes, but patterns of genetic subdivision were not detected in the nuclear markers. Thus, historic isolation of refugial populations has led to prominent mitochondrial lineages but not to cryptic species because populations are interbreeding. The lack of obvious pairwise differentiation among individual populations suggests that D. cephalotes is capable of dispersing between different headwater streams, presumably by flying female imagines. Our data suggest that, within a few generations, D. cephalotes should be capable of reaching and recolonizing restored habitats at a regional scale. Therefore, D. cephalotes is likely to be a valuable and reliable indicator species for monitoring the success of stream restorations if found in restored freshwater sections from which it was absent prior to restoration. Dispersal of D. cephalotes is a valuable case study, but restoration success per se must be estimated based on a wider range of taxa.
We thank Uwe John (AWI Bremerhaven) for carrying out the 454 sequencing. Andrey Rozenberg and Philipp Brand (Ruhr University Bochum) kindly helped with bioinformatic analyses. We also thank the members of the Evo Eco Journal Club and 2 anonymous referees for helpful suggestions that improved this manuscript. This work was supported in part by the GeneStream project ( http://GeneStream.de), which is funded by the Kurt Eberhard Bode Foundation within the Deutsches Stiftungszentrum, and by DFG grant no. HE 2764/2-1.