Genetic diversity in the Alpine flatworm Crenobia alpina

The freshwater flatworm Crenobia alpina (Platyhelminthes, Tricladida, Planariidae) lives almost exclusively in cold springs and crenal streams and possesses only limited dispersal ability. In this study fragments of the COI and 18S rRNA genes were used to estimate genetic divergences among 37 C. alpina populations from the European Alps. Phylogenetic analyses revealed five geographically and genetically distinct groups and at least 10 distinct lineages of C. alpina across the European Alps. Our study suggests that C. alpina represents a complex of numerous cryptic species. Speciation (allopatric and/or sympatric) may have been facilitated by the orogenetic activity of the Alps and the high habitat specificity.


Introduction
Although the Alps are a relatively young mountain range (roughly 80 million years old; Veit, 2002), numerous geological and climatic events have repeatedly influenced the geographical structure of these mountains (Willet et al., 2006), with strong impacts on diversity and the distribution of Alpine life (Thienemann, 1950;Holdhaus, 1954;de Latin, 1967;Schönswetter et al., 2005).In this study we analyzed molecular sequence divergence of the freshwater flatworm Crenobia alpina (Dana, 1766) across the European Alps.C. alpina occurs in most parts of Europe, including Corsica and other Mediterranean islands (Brändle et al., 2007).Outside the Alps, this flatworm lives almost exclusively in cold springs (Durance and Ormerod, 2010;Roca et al., 1992).In the Alps, C. alpina also occurs in crenal streams and at the bottom of Alpine lakes at up to 3000 m a.s.l.(Thienemann, 1950).Springs and crenal streams at high elevations are extreme environments (e.g., constant low water temperature and nutrient availability) that may impose stabilizing selection that results in morphological convergence during speciation (morphological stasis; Rothschild and Mancinelli, 2001;Nevo, 2001;Lefébure et al., 2006a;Bickford et al., 2006).Thus, diversity in extreme environments is often underestimated by classic taxonomic methods (e.g., Trontelj et al., 2009).This is particularly true for taxa with few distinct morphological characters, such as flatworms (e.g., Rader et al., 2017).Our previous genetic study (Brändle et al., 2007) revealed two distinct lineages of C. alpina distributed in the low mountain ranges of Germany and provided first evidence for the strong divergence among populations of C. alpina along the northern margin of the Alps.In the current study our aims were twofold: (i) to characterize the genetic divergence of populations of C. alpina across the European Alps using gene sequences (mitochondrial CO I, 18S rRNA) and (ii) to analyze whether genetic lineages of C. alpina show a spatial structuring across their Alpine range.With our study we contribute towards a better understanding of the role of mountain uplift for evolutionary processes as well as for a better understanding of the role of mountains in structuring the distribution of biodiversity.Such an aim requires a large number of studies on a broad range of taxa with diverging life history traits.By using a flatworm species we provide information from a understudied lineage.

Organism and sampling
In August 2006 and May 2007 10-20 individuals of C. alpina were sampled from 37 springs and headwaters across the Alps (Fig. 1, Table S1 in the Supplement).To gain a first insight into the genetic diversity of C. alpina only one individual from each population was chosen for sequencing.In addition, four individuals from four populations in Corsica,  2).Groups with distinct geographic structuring: southwestern Alps -white diamonds; western Alps -white triangles; northwestern Alps -black triangles; southern Alps -black circles; eastern Alps -black squares.Groups without distinct geographic structuring: mixed clade I -white circles; mixed clade II -white squares.Haplotypes indicated by grey triangles could not be assigned to a geographic group.three individuals from three populations in Sweden, one individual from a British population and one individual of Crenobia alpina montenigrina sampled in Macedonia from a tributary of Lake Ohrid were included in our study (for more details of sampling localities see Table S1 in the Supplement).These individuals from populations outside the Alps were considered to get a first impression of the large scale genetic diversity of C. alpina in Europe.One individual of Polycelis felina sampled in Switzerland was used as outgroup.Specimens were preserved in 99 % ethanol and stored at −23 • C. All individuals analyzed in this study are listed in Table S1.Since C. alpina is mostly restricted to springs and headwaters each sampling location was defined as a separate population.

DNA isolation and amplification
The DNeasy ™ Tissue Kit (Qiagen) was used to extract DNA following the protocol provided for animal tissues.For the mitochondrial cytochrome oxidase subunit I (CO I) gene the primers Pr-a2 (5'-AGC TGC AGT TTT GGT TTT TTG GAC ATC CTG AGG T-3'; Bessho et al., 1992) and CO13B (5'-AAG TGT TGN GGR AAR AAN GT-3'; Telford et al., 2000) were used for amplification.These primers increased the length of the sequenced fragment considered in our previous study (see Brändle et al., 2007) from 382 to 555 bp.18S rRNA gene fragments with a length of 544 bp were amplified using the primers 18S1 (5'-TAC TGT TGA TCC TGC CAG TA-3'; Kuznedelov et al., 1996) and 18S2 (5'-ATT ACC GCG GCT GCT GGC ACC-3'; Kuznedelov et al., 1996).The fragments were amplified in a thermocycler (Eppendorf Mastercycler, Hamburg, Germany).The PCR conditions for the COI gene consisted of 5 min initial denaturation at 95 • C, followed by 31 cycles of 30 s denaturation at 94 • C, 50 s annealing at 50 • C and 1 min extension at 72 • C, followed by 10 min final extension at 72 • C. The PCR conditions for the 18S rRNA gene consisted of 5 min initial denaturation at 94 • C, followed by 31 cycles of 1 min denaturation at 94 • C, 70 s annealing at 55 and 2.5 min extension at 74 • C, followed by 10 min final extension at 74 • C. Products were purified using the MinElute ™ PCR purification kit (Qiagen) following the protocol provided.The purified PCR fragments were sequenced in both directions by a commercial company (GENterprise GmbH, Mainz, Germany).All haplotypes were submitted to GenBank (accession number KY569117-KY569175; see Table S1).

Phylogenetic analyses
Partial COI gene and 18S rRNA gene sequences were aligned using BioEdit 7.0.5.3 (Hall, 1999).In a number of sequences particularly within the COI fragment, some positions showed ambiguities and were discarded.These exclusions reduced the length of the COI gene fragment to 450 bp and those of the 18S RNA gene fragment to 531 bp.Thus, the total length of the combined fragments used for phylogenetic inference was 981 bp.We also translated the COI gene sequences into amino acid sequences according to the flatworm mitochondrial genetic code (Nei and Gojobori, 1986) using MEGA 5.05 (Tamura et al., 2011).Except for one position in the sequence of C. alpina montenigrina that was translated into an internal stop codon, no insertions, deletions or stop codons were found.In the 18S RNA gaps were coded as missing data.
The number of substitutions and genetic distances between haplotypes of both fragments and the combined data sets were calculated with MEGA 5.05.Spatial structuring of the sampled populations was analyzed by means of correlation analysis of pairwise genetic versus pairwise geographic distances among the sample sites from the Alps (all 37 sampling sites) using patristic distances inferred from the Bayesian tree of the COI gene fragment (see below).Patristic distances were calculated with the program PATRISTIC (Fourment and Gibbs, 2006).We tested this matrix correlation using a Mantel test as implemented in the "vegan" package (Oksanen et al., 2013) in the software R (www.r-project.org).
Phylogenetic relationships among haplotypes were analyzed with MRBAYES version 3.1.2(Ronquist and Huelsenbeck, 2003).To select an appropriate evolutionary model of nucleotide substitution for each nucleotide fragment, the approach of Minin et al. (2003) as implemented in the Perl script DT-ModSel was used.In contrast to the likelihood-ratio tests implemented in MrModeltest 2.2 (Nylander, 2004), this approach of performance-based model selection uses the Bayesian information criterion and considers the relative branch-length error as a performance measure in a consistent framework with penalty for overfitting based on decision theory (Minin et al., 2003).MRBAYES was run for 10 6 generations, submitting the two gene regions separately as well as both together.Since different models of sequence evolution for the different genes were detected (see results), partitions for the two genes have been specified.Trees were sampled every 1000th generation using a random tree as a starting point.Inspection of the standard deviation of split frequencies after the final run indicated convergence of Markov chains (< 0.03).In all analyses, two parallel Markov chain Monte Carlo simulations with four chains (one cold and three heated) were run.The first 10 samples of each run were discarded as burn-in.From the sampled trees, consensus trees were produced using the ALLCOMPAT command in MRBAYES.

Results
For the CO I gene, a total of 38 haplotypes were found among the 46 sequenced Crenobia specimens; 187 sites (33.7 %) were variable and 156 sites were parsimony informative (28.1 %).Of the variable sites, 84.5 % were found at the third codon position, 13.4 % at the first codon position and 2.1 % at the second codon position.The average number of base differences between sequences was 60.3 (±4.49SD).For the 18S rRNA gene, 19 haplotypes were found across the 46 sequenced specimens; 27 (5 %) positions were variable and 17 (3 %) were parsimony informative.The average number of base differences between sequence pairs was 7.21 (±1.45 SD).The combined data set was characterized by 42 haplotypes with 214 variable sites (19.5 %) and 177 (16.1 %) parsimony-informative sites.The average number of base differences per sequence pairs was 69.9 (± 4.59 SD).DT-ModSel selected the HKY+I+G model for the COI gene and the HKY+I model for the 18S rRNA gene (computed on the haplotype files, outgroup excluded).
Phylogenetic analysis of the combined data set revealed distinct geographic groups generally supported by high posterior probabilities (Fig. 2).We named these groups according to their geographical origin.Distinct groups were found in the western, southern and eastern parts of the Alps, whereas no specific group was located in the central parts of the Alps.Populations analyzed from outside the Alps were both strongly different (e.g., Corsica) but also very similar (e.g., Great Britain) from the populations of the Alps.Bayesian trees calculated across the individual genes (Figs.S3 and S4 in the Supplement) revealed a very similar grouping of at least the basic lineages when compared with the Bayesian tree of the combined data set (Fig. 2).The relationships among these lineages, however, differed considerably between the data sets.All basal splits were only weakly supported by posterior probabilities.In an attempt to increase the support for the basal branches, we included the secondary structures of the 18S rRNA gene (for methods, see Additional Methods S5).However, also with this approach, the basal splits were only weakly supported.Although we sequenced only one individual from each population, we found a positive correlation between geographic distances (log 10transformed) and genetic distance among the samples from the Alps (Fig. 3, Mantel test: r = 0.31, P < 0.001, 999 permutations).However, this relationship is low and characterized by considerable scatter.

Discussion
The main results of our study can be summarized as follows: we found (i) considerable genetic diversity within C. alpina in the European Alps with at least 10 genetic lineages arguably cryptic species and (ii) a significant but low spatial structuring of haplotypes/lineages.

Genetic diversity and cryptic species
The distinct lineages recognized in the present study may represent different species.In general, species are considered cryptic if they have been erroneously classified as a single species because of a lack of distinct morphological features.DNA barcoding has revealed cryptic diversity in a variety of taxa (Nieberding et al., 2005;Bickford et al., 2006;Witt et al., 2006;Lefébure et al., 2007;Pauls et al., 2010).Lefébure et al. (2006b) showed that clades of crustaceans diverging by more than 0.16 substitutions per site of the CO I gene as measured by patristic distances are likely to present distinct species.Based on the threshold distance suggested by Lefébure et al. (2006b) a minimum of 10 distinct species of C. alpina may be present across the Alps.Further evidence for putative species was found by analyzing uncorrected pairwise sequence divergences.Uncorrected p-distances of the CO I gene ranged from 0.1 to 14.6 % (±0.03 SD), with an average of 11 %.In a study of the cryptic diversity of Tetramorium ants in the western Palearctic, Schlick-Steiner et al. (2006) found interspecific divergences for the CO I gene ranging from 4.0 to 10.6 %.Liu et al. (2003) investigated cryptic species diversity within the western American springsnail Pyrgulopsis micrococcus and found uncorrected p-distances of the CO I gene among lineages from 3.69 to 10.6 %.Our phylogenetic analyses provided some evidence for distinct species also within C. alpina.Descriptions of new species in flatworms, however, will require studying the internal anatomy in particular the copulatory apparatus of planarians.But note that body form, reproduction mode and morphology of the reproduction apparatus of the Turbellaria vary considerably depending on environmental conditions and food supply (e.g., Steinböck, 1942).As a consequence, classical taxonomy may have overlooked species in the past leading to an underestimation of the diversity of lineages (see also Rader et al., 2017).In that regard we hope that genetic characterization could guide future morphological studies.Pfenninger and Schwenk (2007) showed that cryptic species are evenly distributed among the major metazoan taxa.This suggests that a simple body plan is not a suf-  S1 and S2).Posterior probabilities (> 0.8) are shown on the nodes.Haplotypes were named according to their geographical location, and the names are followed by the country code in parentheses (Supplement Table S1).Genetically and geographically distinct groups are indicated by braces referring to the geographic position within the Alps on the right side of the figure.For details about sampling localities, see Supplement.
ficient explanation for an increased occurrence of cryptic species within a taxonomic group.Some authors (e.g., Bickford, 2006) have suggested a link between extreme environmental conditions and cryptic species.Extreme environments impose stabilizing selection on changes in morphology during speciation, resulting in morphological stasis (Rothschild and Mancinelli, 2001;Nevo, 2001;Lefébure et al., 2006a;Bickford, 2006;Lefébure et al., 2007).C. alpina is in general restricted to extreme habitat conditions of cold springs and crenal streams.Harsh conditions occurring in such environments and emerging biological properties, such as general inability to tolerate temperatures above 13-15 • C and preference of temperatures slightly above 0 • C (Steinböck, 1942;Reynoldson, 1953;Beauchamp and Ullyott, 1932;Köhler, 1932, in Wright, 1974), might limit morphological divergence.Evidence for morphological stasis of groundwater amphipods (Lefébure et al., 2007;Trontjl et al., 2009) and spring-inhabiting snails (Benke et al., 2009) has been found.The concept that extreme environmental conditions favor the evolution of cryptic species still requires further evaluation.Pfenninger and Schwenk (2007) have already shown that there are no differences in the occurrence of cryptic species among different biogeographical regions, which raises doubt about the importance of environmental conditions for cryptic diversity.

Spatial structuring
Given that we propose the existence of cryptic species and since our sampling scheme is based on one individual per site only, we cannot detect multiple species at one location.
In fact, it is surprising that we still found moderate spatial structuring in our data.According to our phylogenetic analysis of C. alpina, roughly five biogeographic groups could be distinguished: (i) western Alps I, (ii) southwestern Alps, (iii) eastern Alps (iv) western Alps II and (v) southern Alps (Fig. 1).In the central parts of the Alps, we found no clear geographical structuring of C. alpina.The detected spatial structuring fits with previous findings for plants (Schönswetter et al., 2005) and animals (e.g., Pauls et al., 2006) and thus might be a promising hint for future studies with more samples per site and additional markers.With such approaches it would be possible to retrieve a more complete and convincing picture of the inter-and intraspecific divergence and to test whether closely related cryptic species show vicariant patterns as documented between flatworm genera (Wright, 1974;Snook and Milner, 2001;Roca et al., 2006).Furthermore, due to the strong environmental link to cold springs it is expected that little gene flow occurs between flatworm populations (Brändle et al., 2007).The moderate geographic structuring revealed in our study may support this finding.However, the well-known potential for asexual reproduction in C. alpina (Reynoldson and Young, 2002) and founder ef-fects may also have a strong influence on the geographic structuring.
Our study provides some evidence that the orogenic history of the alps and presumably the climatic variability during the Quaternary promoted genetic diversity of C. alpina.Although our study represents only a first insight into the evolutionary history of C. alpina in the Alps, we demonstrated that it is well worth having a closer look at species that apparently show less variability in terms of morphology, habitat choice and ecology (Rader et al., 2017).Judging whether the detected genetic lineages represent different species will be a challenge of future studies.

Figure 1 .
Figure 1.Map of the Alps with locations of sampled Crenobia alpina populations.Different symbols refer to the grouping revealed by the phylogenetic analyses (see Fig. 2).Groups with distinct geographic structuring: southwestern Alps -white diamonds; western Alps -white triangles; northwestern Alps -black triangles; southern Alps -black circles; eastern Alps -black squares.Groups without distinct geographic structuring: mixed clade I -white circles; mixed clade II -white squares.Haplotypes indicated by grey triangles could not be assigned to a geographic group.

Figure 2 .
Figure 2. Bayesian phylogeny of Crenobia alpina haplotypes based on mitochondrial COI and nuclear 18S rRNA gene fragments (overall 981 bp).Tip labels refer to locations sampled and analyzed (Supplement TablesS1 and S2).Posterior probabilities (> 0.8) are shown on the nodes.Haplotypes were named according to their geographical location, and the names are followed by the country code in parentheses (Supplement TableS1).Genetically and geographically distinct groups are indicated by braces referring to the geographic position within the Alps on the right side of the figure.For details about sampling localities, see Supplement.

Figure 3 .
Figure3.Relationship between pairwise genetic distances of the COI gene (patristic distances from tree shown in Fig.S3) and pairwise geographic distances of Crenobia alpina populations across the European Alps.Note that the abscissa (geographic distance) is log 10 transformed.