2 Defining units for conservation management for the European tree frog
2.6 Acknowledgement
This research was supported by grants from the German Federal Environmental Foundation (DBU), Heidehof-Stiftung, and “Hans-Schiemenz-Fonds“ - Deutsche Gesellschaft für
Herpetologie und Terrarienkunde (DGHT). I thank the following Nature conservation authorities for permission for tree frog collection: the biosphere reserve Niedersächsische Elbtalaue, Kreis Minden Lübbecke, Kreis Steinfurt, Landkreis Diepholz, Landkreis Gifhorn, Landkreis Lüneburg, Landkreis Osnabrück, Landkreis Stade, Landkreis Uelzen, Region Hannover, Stadt Wolfsburg, Sachsen Anhalt. I am especially grateful to Annika Ruprecht,
Christina Akman, Frank Weihmann, Günter Krug, Heike Pröhl, Irena Czycholl, Ivonne Meuche, Jana Kirchhoff, Kim Jochum, Matei Balborea, Michael Weinert, and Wiebke Feindt for help during field work. Finally I thank our technician Sabine Sippel for her assistance in the molecular lab.
3
Phylogeographic structure of the European tree frog
(Hyla arborea) in its German distribution area
3.1 Abstract
Knowledge about the existence of different genetic lineages within endangered species is important for conservation management. To assess the phylogeographic structure of the European tree frog across its distribution area in Germany, 372 individuals were sampled at 31 sites and sequence analyses of a mitochondrial gene fragment (cytochrome b) and analyses of eight microsatellite loci were carried out. Sequence divergence between sample sites was low varying between 0 and 0.4 % (overall: 0.2 %) and no distinct genetic lineages were found. Nonetheless, a clear North-South partitioning could be revealed by both molecular markers with the Central German Uplands as probable barrier. Furthermore, the influence of the major rivers such as Elbe, Rhine, and Danube on the phylogeographic structure could be revealed. Concerning future conservation measures, the identified genetic structures should be considered, especially for the choice of individuals if resettlements should be necessary.
3.2 Introduction
The current genetic structure of many species was fundamentally shaped by the recent glacial periods. Depending on the location and number of glacial refugia and the route of
recolonisation of the continent, distinct genetic lineages of a species could have evolved.
Especially for threatened species for which conservation measures are planned or are already implemented, the boundaries of different genetic lineages and potential contact zones are important to recognize to decide about strategies to avoid accidentally mixing distinct clades.
Besides the general aim to maintain the species genetic diversity and historic integrity, the fitness of hybrid offspring can suffer reduced fitness due to endogenous selection by melding distinctly co-adapted gene complexes (Harrison, 1993).
Phylogeographic analyses, mostly based on sequencing of mitochondrial DNA (cyt b, d-loop, COI), can reveal such lineages. In Europe different genetic lineages were described for several species (Taberlet et al., 1998 e.g. Ursus arctos, Crocidura suaveolens, Chorthippus parallelus). Also for amphibian species distinct genetic lineages were described. For example two distinct Rana temporaria lineages were described in Eastern and Western parts of Europe forming a contact zone in northern Germany. Furthermore, evidence was found for an Irish glacial refugium (Palo et al., 2004; Schmeller et al., 2008; Teacher et al., 2009). For the fire salamander (Salamandra salamandra) an Eastern and Western postglacial lineage were described (Steinfartz et al., 2000) forming contact zones in Western Germany (Weitere et al., 2004)
The presence of such potential distinct lineages or contact zones needs to be examined for the European tree frog, Hyla arborea, a species that showed long-term decline in much of its Western European distribution, mainly caused by habitat loss, fragmentation and
degradation. Conservation measures such as habitat restoration and population resettlements are conducted. Stöck et al. (2008) made a first attempt with a Europe-wide molecular
analysis. Since their purpose was a wide-ranging phylogenetic study, only few sample sites in Germany were included (four individuals from three sites in Germany) and therefore
important genetic structures could have been overlooked. Therefore, my interest is, with a more extensive sampling network, to analyse the phylogeographic structure of tree frogs in Germany. My main question is if there are two or more genetic lineages of the European tree frog in different areas in Germany. The results will be discussed in the context of already
published mitochondrial data of other European sample sites and the importance of the results for conservation management.
3.3 Material and methods
3.3.1 Sample collection and preparation
Thirty-one sites were sampled spanning a coarse net across the tree frog distribution in Germany. In total 372 individuals were sampled with 1 - 22 (mean: 12) individuals per sample site (see Table 3.1). Because of a preliminary study on the relationship of European tree frogs in Lower Saxony the sampling is more detailed in this area. Genetic material was collected by tips of tadpole tails and by buccal swabs of adult frogs. The adults (mostly males) were collected from the choruses during the breeding season in spring 2005, 2008 and 2009. Tadpoles were sampled in summer 2007. DNA from the tail clips was fixed in 99 % ethanol and extracted using a proteinase K digestion followed by a Phenol-Chlorophorm protocol (Sambrook et al., 1989) and stored at -20 °C. From the buccal swabs DNA was extracted with an Invisorb Spin Swab Kit (Invitek) following the manufacturer’s protocol and stored at -20 °C.
Fragments of 900 bp of cytochrome b (cyt b) of 1 - 20 individuals from each sample site, except KZ and BH in the Hannover population, were amplified via PCR using the primers MVZ 15-L (5′- GAACTAATGGCCCACACWWTACGNAA -3′) and Cytb AR-H (TAWAAGGGTCTT CTACTGGTTG) from Moritz (1992) and Goebel (1999). The PCR reaction (25 µl) consisted of 20 - 100 ng DNA, 1 µl of each primer (10 µM), 0.8 µl dNTP’s (10mM 5PRIME), 2.5 µl 10x advanced Buffer (5PRIME), 1.25 U Taq DNA Polymerase (5PRIME), and 17.45 µl H2O. PCR conditions were as follows: an initial denaturation at 94 °C for 3:00 min; 35 cycles at 94 °C for 45 s, annealing temperature of 50 °C for 45 s, extension at 65 °C for 1:00 min. The PCR products were sent to the Macrogen Company (Seoul, South Korea) for purification and sequencing with an ABI3730XL genetic analyzer (Applied Biosystems).
Table 3.1: Overview of sample sites. a: samples from adult frogs, t: samples from tadpoles, Ho: observed heterozygosity, He: expected heterozygosity, SD: standard deviation, NA: mean number of alleles over all loci, h: haplotype diversity, π: nucleotide diversity, N: number of sampled individuals for microsatellite respectively sequence analyses
Since preliminary analyses revealed low variation of cyt b sequences I included the analysis of eight species specific microsatellite loci (WHA1-9, WHA1-20, WHA1-25, WHA1-60, WHA1-67, WHA1-103, WHA1-104, and WHA1-140) previously isolated by Arens et al. (2000). The microsatellites were amplified for 1-22 individuals from each sample site following the authors’ protocol, except for the annealing temperature for WHA1-20, which was changed to 64.6 °C. PCR products were genotyped using the capillary sequencer MegaBace 1000 (Amersham Bioscience) and ABI 3500 (Applied Biosystems). Allele scoring was performed using the corresponding software Genetic Profiler v2.2 and GeneMapper v4.1.
3.3.2 Statistical analysis
3.3.2.1 Analysis of mtDNA in Germany
Both directions of the cyt b sequences were assembled using the computer software SeqMan™ II (DNASTAR, Inc., Konstanz, Germany). Multiple sequence alignments were performed in MEGA 5 (Tamura et al., 2011) using the Muscle algorithm (Edgar, 2004) and all variable sites were confirmed by visual inspection of the chromatograms. The same program was used to calculate p-distances between sample sites (Tamura et al., 2004). A haplotype network of the cyt b data set was constructed via the statistical parsimony analysis of the program TCS 1.21 (Clement et al., 2000) using the default settings. Haplotype diversity (h) and nucleotide diversity (π) (Nei, 1987) were determined with ARLEQUIN v. 3.11
(Excoffier et al., 2005).
3.3.2.2 Analysis of microsatellites in Germany
Microsatellite-data were checked for null alleles, stuttering and allelic dropout using MICRO -CHECKER (Van Oosterhout et al., 2004). The program FSTAT v. 2.9.3 (Goudet, 1995) was used to test for genotypic disequilibrium of all pairs of loci in each sample. For each sample site and locus the observed and expected heterozygosity (Nei, 1987) and deviation from Hardy-Weinberg equilibrium (HWE) (Guo and Thompson, 1992) were determined with ARLEQUIN
v. 3.11 (Excoffier et al., 2005). GENEPOP v. 4.1 (Rousset, 2008) was used to test for a global deviation from HWE in each sample site. Genetic differentiation between the sample sites
was calculated as pairwise Dest values (Jost, 2008) using the R package DEMEtics (Gerlach et al., 2010). Significance was calculated by 10,000 bootstraps.
A Bayesian clustering model performed with the program STRUCTURE 2.3.3 (Pritchard et al., 2000) was used to infer genetic clusters. The aim of this method is to define clusters of individuals on the basis of their genotypes at multiple loci using a Bayesian procedure. It attempts to find population clusters by reducing linkage disequilibrium and deviations from the Hardy-Weinberg equilibrium within inferred clusters. The user specifies a priori the number of population clusters (K) and estimates the log likelihood Pr(X|K) for this model. For finding the most likely number of genetic clusters the log likelihood Pr(X|K) is always
calculated for a series of K values.
All STRUCTURE runs used 200,000 iterations after a burn-in period of 50,000. Because of the large distances between most sample sites, I used the assumption of the no-admixture model and independent allele frequencies. Fifteen runs were performed for each K. The range of possible Ks tested was from 1 to 31, according to the number of sampled breeding sites. I calculated the average log likelihood Pr(X|K) (given by the estimated Ln Prob of data = Ln P(D) in the software result output) for each K across all runs. Since it is not always
straightforward to detect the true number of K, I included the ΔK statistics proposed by Evanno (2005) using Structure Harvester v.0.56.4 (Earl and vonHoldt, 2012). Furthermore, STRUCTURE provided values of allele-frequency divergence among revealed clusters (net nucleotide distance), computed by using point estimates of P.
Because of low sample size populations BLO, AGM, HAK, HKT and SAZ have been excluded from most microsatellite analyses. Sample sites SB and QU were included,
nonetheless results should be regarded with caution. To assess significances in this study I applied sequential Bonferroni corrections (Rice, 1989) to all multiple comparisons.
3.3.2.3 Analysis of mtDNA in the European context
To put the phylogeographic structure found in Germany in a species-wide context, I compared my data with further 90 cyt b sequences of 30 sample sites of Hyla arborea in Europe which were already published by Stöck et al. (2008; 2011). Another haplotype network was built on the basis of 851 bp fragments.
3.4 Results
3.4.1 Analysis of mtDNA in Germany
I revealed 26 haplotypes of the cytochrome b fragment (Appendix 1) which differed by 28 variable sites and 23 parsimony informative sites (Figure 3.1; for detailed information on individual-haplotype-assignment see Appendix 2). Most haplotypes were unique to one sample site except haplotype Hy-1, Hy-2, Hy-5, Hy-12, and Hy-20. Hy-5 (red) was the main haplotype showing a broad distribution almost over the complete sampling area. All other haplotypes differed by only single mutation steps. The second most common haplotype Hy-1 (blue) and its derivates Hy-6, Hy-9, Hy-21, and Hy-22 were restricted to the Northern part of Germany. Hy-2 (green) is shared by five sample sites near the river Elbe and Hy-12 (yellow) by two sample sites near the river Danube. The two sample sites sharing the haplotype Hy-20 (violet) are both located in the Middle East of Germany. (Figure 3.2)
P-distances among sample sites were low, varying between 0 and 0.4 % (Appendix 6).
The overall distance was 0.2 %. The highest estimates of mtDNA diversity were found in SAZ (h = 0.83, π = 0.11 %) and OL (h = 0.79, π = 0.21 %). The lowest values were found in BLO, HAK, QU, EK, KZ, and AGM (all: h = 0, π = 0 %).
Figure 3.1: Haplotype network of 26 distinct haplotypes of cytochrome b of Hyla arborea (900 bp) in Germany.
Each haplotype is represented by one circle. The size of the circles corresponds to the haplotype frequency.
Lines between haplotypes denote mutational steps between sequences; black nodes denote inferred intermediate haplotypes between observed haplotypes. AN, BA etc. are abbreviations sample sites where the haplotype was found. Haplotypes shared by different sample sites are marked in strong colours; haplotypes which are present in one sample site only are left blank. For pattern illustration haplotypes with light blue framing mark haplotypes originating from the blue haplotype
Figure 3.2: Haplotype distribution. Physical map of Germany with distribution of cyt b haplotypes in the sampling area in Germany. Each haplotype is represented by one colour corresponding to the colours in the haplotype network (Figure 3.1).Green points denote the current distribution of Hyla arborea on the basis of TK25 modified from: Report on the main results of the surveillance under article 11 for annex II, IV and V species (Annex B) (http://cdr.eionet.europa.eu).
3.4.2 Analysis of microsatellites in Germany
The eight microsatellite markers examined were polymorphic with eight to twenty-two alleles per locus. The analysis using Micro-Checker uncovered signs of null alleles for the locus WHA1-67 in the sample site KO, for the locus WHA1-104 in the sample site KAS and for the locus WHA1-140 in the sample sites SW and OVH. As null alleles for the three loci were found only at single sample sites, I did not adjust for null alleles. Furthermore this analysis revealed no evidence for large allele dropout or scoring errors due to stuttering.
Deviation from Hardy-Weinberg-Equilibrium was found for WHA1-60 with a significant heterozygosity excess in the sample sites KH, BH and AN. For WHA1-104 a deficiency was found in KH and KAS. The global test for HWE over all loci in each population resulted in no significant deviation. No Linkage (genetic) disequilibrium was found between any pair of loci.
Since Berset-Brändli et al. (2007) found the locus WHA1-60 to be sex linked with a suppressed recombination in males, I tested its influence on the outcome of all analyses.
Expected heterozygosity values did not change remarkably after excluding WHA1-60 (mean He locus WHA1-60 included: 0.73; excluded: 0.72). An influence of WHA1-60 on the results of all other analyses was not evident. Therefore, I decided to keep this locus in the analyses.
The highest genetic diversity, expressed as expected heterozygosity, was found in KAS (He = 0.83) and KOB (He = 0.82). The lowest values were found in WK (He = 0.60) and MOZ (He = 0.63).
Pairwise genetic distances measured as Dest values were found to vary between 0 and 0.7 (Appendix 6). Whereas, except the two closest sites KZ and KO, seven other comparisons were found to be not significantly different. Since all these comparisons include SB – a site with only five sampled individuals – these values should be regarded with caution. In general the highest Dest values (> 0.6) were found when comparing the southern sample sites RAV, BGH, and ISM with the Northern sample sites.
For the Bayesian analysis a first peak in ΔK is found at K = 2 (Figure 3.3) separating the southern sample sites KOB, FS, KAS, OFF, RAV, BGH, ISM and additionally MOZ from the Northern sites. Since MOZ shows in the further separation (K = 4) a higher relationship to the northern clusters, as indicated by lower net nucleotide distances (Appendix 8), I
conducted no hierarchical analysis.
A sign for a first plateau in Ln P(D) is found at K = 4. One of the four genetic clusters (red) is found in the South of Germany. Two Clusters (blue and green) were found in the North of Germany. The blue cluster encompasses large parts of the Northwest and the green cluster is mainly distributed in the Northeast. A fourth sharply separated cluster is found for the sample site MOZ (grey) (Figure 3.4).
Further subclustering is indicated by another peak at K = 8 in ΔK method and larger Ln P(D) values. The southern genetic cluster is further separated into a yellow cluster along the Danube River, an orange cluster in the Upper Rhine Plain, and a red cluster along the Main River. In the North further separation is revealed by a cluster occurring mainly in the sample sites in the West of Hannover (dark blue) and a cluster in light blue occurring mainly in the North eastern part of Lower Saxony. In light green, a cluster separated from the former, mainly occurring in SW and OL in the Northeast (detailed barplots are given in Appendix 7).
Again MOZ is sharply separated from all other sites by the grey cluster. This is elucidated furthermore by relatively high network distances in comparison to the other clusters (Appendix 8 and 9).
Figure 3.3:Mean values of estimated Ln probability of data (LnPD) for each K (a) and delta K (b)
Figure 3.4: Distribution of distinct genetic clusters revealed by Structure analysis of microsatellites for K = 4 (left) andK = 8 (right). Each cluster is represented by a different colour. Green points give the distribution of Hyla arborea
3.4.3 Analysis of mtDNA in the European context
Including 90 further cyt b sequences from other sample sites in Europe (Stöck et al., 2008;
Stöck et al., 2011) in the network analysis, I found 48 haplotypes which are still very similar (Figure 3.5). Some variation is found in Greece and on Crete. Nonetheless, haplotypes of this branch differ only in few base pairs.
Haplotype Hy-5 (red) is widely distributed from Western France to Albania and Romania. Hy-1 (blue) is found in the Netherlands extending the structure found in North Germany further into the Northwest (Figure 3.6).
Figure 3.5: Haplotype network of 48 distinct haplotypes of cytochrome b of Hyla arborea (851 bp) in Europe.
Each haplotype is represented by one circle. The size of the circles corresponds to the haplotype frequency.
Lines between haplotypes denote mutational steps between sequences; black nodes denote inferred intermediate haplotypes between observed haplotypes. BE, DE, GR etc. denote code of country where the haplotype was found. Haplotypes shared by different sample sites are marked in strong colours; haplotypes which are present in one sample site only are left blank. For pattern illustration haplotypes with light blue framing mark haplotypes originating from the blue haplotype and haplotypes with light brown framing mark the divers branch found in Greece.
Figure 3.6: Distribution of cyt b haplotypes in Europe. Each haplotype is represented by one colour corresponding to the colours in the haplotype network (Figure 3.5). Close sites in Belgium, Croatia, the Netherlands, Romania, and on Crete are grouped together.
3.5 Discussion
I conducted a large scale molecular analysis on the European tree frog in its German
distribution area to assess the potential for the presence of distinct genetic lineages thought to have evolved due to the postglacial recolonisation of the continent. Sequence divergence between sample sites was low. Nonetheless, a clear North-South partitioning could be
revealed by both molecular markers. Within the phylogeographic structure, the influence of higher mountain ranges and the major rivers becomes apparent.
3.5.1 Distinct genetic lineages in the European tree frog?
For several European species distinct genetic lineages could be detected. They evolved by the separation in different glacial refugia and the subsequent recolonisation of the continent. In amphibians for Rana temporaria an Eastern and Western lineage was found (Palo et al., 2004). The mean interlineage divergence between the cyt b haplotypes in the two clades was 3.2 %. For Rana arvalis two main clades were detected differing by 3.6 % cyt b sequence divergence (Babik et al., 2004). It was inferred that these clades survived several glacial cycles. One of these clades showed two further subclades which arose presumably during the last glaciation. They differed by 1 % sequence divergence.
For the European tree frog however such a strong genetic structure could not be detected. Cyt b sequence divergence (p-distances) between sample sites in this study was low varying between 0 and 0.4 % (overall: 0.2 %). Indicating that no different genetic lineages are apparent neither in the German, nor in the European distribution area of the European tree frog. So far this supports the hypothesis of Stöck et al. (2008). Since they found European wide homogeneity of mtDNA but higher diversity of nuDNA in the Balkan region they supposed the spread of a single mtDNA lineage from a potential Pleistocene refugium in the Balkan region. The star like haplotype structure found in this study suggests a rapid
postglacial colonization of the continent. The results of an mtDNA mismatch distribution analysis of Stöck et al. (in press) on the smaller sample set also pointed to a recent and rapid expansion in this species.
3.5.2 Phylogrographic structures of the tree frog in Germany
Although haplotypes are very similar, and mostly diverged in only one base pair, some phylogeographic structure could be detected in the German distribution area. For example a clear North-South partitioning could be revealed by both molecular markers. Haplotype Hy-1 (blue) is widely distributed in Northern Germany and adjacent areas in the Netherlands.
Present only in their own sample site, but originating from 1, are the four haplotypes Hy-6, Hy-9, Hy-21, Hy-22 (pastel blue Figure 3.1). Interestingly, the border of the distribution of Hy-1 and its descendants (blue branch) is congruent to the border of the Central German Uplands (Figure 3.2) which seemed to have had a major role as migration barrier. It is likely that the mutation between Hy-5 and Hy-1 happened during the colonisation into the North. A rapid distribution in the almost barrier free North German Plain could have led to the
Present only in their own sample site, but originating from 1, are the four haplotypes Hy-6, Hy-9, Hy-21, Hy-22 (pastel blue Figure 3.1). Interestingly, the border of the distribution of Hy-1 and its descendants (blue branch) is congruent to the border of the Central German Uplands (Figure 3.2) which seemed to have had a major role as migration barrier. It is likely that the mutation between Hy-5 and Hy-1 happened during the colonisation into the North. A rapid distribution in the almost barrier free North German Plain could have led to the