1 General introduction
1.4 Aims of the study
1.4.2 Management Units in Lower Saxony and adjacent areas
In Lower Saxony the current distribution of the European tree frog is very patchy with some larger occurrences in the lowlands. Severe declines have been observed mainly in the second part of the last century (Manzke and Podloucky, 1995). In some places measures for
conservation management have already been implemented and initial success has been observed (e.g. Clausnitzer, 2004 (Celle); Buschmann et al., 2006 (Steinhuder Meer); LaReG, 2007 (Braunschweig); Richter and Mügge, 2012 (Diepholz)). Units for conservation
management need to be delineated for the European tree frog in order to support these conservation activities.
The second aim of this study is to perform a large scale conservation genetic survey of the European tree frog across its distribution in Lower Saxony and adjacent areas and to obtain insight into contemporary as well as historical processes. My special interest is to assess genetic diversity and to define units for conservation management.
2
Defining units for conservation management for the
European tree frog (Hyla arborea) in Lower Saxony
and adjacent areas
2.1 Abstract
The European tree frog Hyla arborea has suffered from dramatic population declines in the last decades and has therefore been categorised as threatened in many Red Data lists. In Lower Saxony in Germany the current distribution of the tree frog is very patchy with some main occurrences in the lowlands. For supporting effective conservation measures this study aims to assess genetic diversity and to define units for conservation management.
Across the tree frog distribution area in Lower Saxony and adjacent areas 237 individuals were sampled at 14 sites (~ 3 - 250 km apart from each other). All samples were genotyped with eight microsatellite loci and twelve sites were sequenced for an mtDNA cytochrome b fragment.
While all but one of the microsatellite pairwise Dest and FST values showed significant genetic differentiation (Dest: 0 - 0.46, FST: 0 - 0.18), Bayesian analyses indicated common structures forming seven distinct genetic clusters. The cytochrome b haplotype distribution highlights the former connection of the currently fragmented populations along the river Elbe.
Since genetic diversity was relatively high, each of the sampled tree frog occurrences should have the potential to recover to a stable population size when applying appropriate local conservation measures. For new resettlement projects identified genetic structures should be considered for the choice of source populations. Where possible, it would be preferable to reconnect originally linked occurrences that are now separated in different conservation unites due to habitat fragmentation and genetic drift.
2.2 Introduction
Amphibian populations around the world are seriously affected by severe declines in the last decades. In Europe habitat loss, fragmentation and degradation are the most significant threats to amphibians (Temple and Cox, 2009). For most of these species intensive conservation efforts are needed to prevent them from further decline and to regain a favourable
conservation status.
Stabilization of weakened populations can be achieved e.g. by improving the habitat, constructing new breeding ponds and connecting populations to stable networks (e.g. Hyla arborea: Tester and Flory, 2004; Bombina bombina: Brockmüller and Drews, 2009; Triturus cristatus and Pelobates fuscus: Rannap et al., 2009). However, in cases where extreme fragmentation isolates populations, reconnection is difficult. Highly isolated populations already suffering from inbreeding and low genetic diversity can be strengthened by
introducing translocated individuals from other populations with the aim of increasing fitness and genetic diversity (e.g. Vipera berus: Madsen et al., 1999). This is a difficult task to undertake because mixing different gene pools can in the best case result in a higher fitness of the offspring (hybrid vigor), but it could also result in outbreeding-effects with fitness
depression (hybrid breakdown) in subsequent generations which may drive the population into further decline. For example the introduction of individuals into a small inbred
population of Florida panthers (Puma concolor coryi) led to a higher survival rate of hybrid offspring and helped to recover the population (Pimm et al., 2006). On the other hand, mixed-source reintroductions of slimy sculpins (Cottus cognatus) have led to outbreeding depression in second-generation descendents. In this case, source populations were genetically
differentiated by an FST of 0.32 (Huff et al., 2011).
Therefore, to minimize negative effects when translocations are necessary or mixed-source introductions are planned, it is essential to reveal genetic structures for effective
species conservation management. Potential hidden barriers and units for conservation need to be delineated.
Units for conservation management have been defined by using genetic analysis for several endangered species such as koalas (Lee et al., 2010), harbour porpoises (Wiemann et al., 2010) or Larch Mountain salamander (Wagner et al., 2005). These studies showed that in all these species limited migration due to natural and anthropogenic barriers formed genetic
distinct units. For purposes of conservation management the studies recommended all handling each unit individually.
The European tree frog is a species that showed long-term decline in much of its Western European distribution, mainly caused by habitat loss, fragmentation and degradation.
In Lower Saxony in Germany the current distribution of the tree frog is very patchy with some main occurrences in the lowlands (Figure 1). Severe declines have been observed mainly in the second part of the last century (Manzke and Podloucky, 1995). At some places measures for conservation management are already implemented and first successes have become apparent (e.g. Clausnitzer, 2004 (Celle); Buschmann et al., 2006 (Steinhuder Meer);
LaReG, 2007 (Braunschweig); Richter and Mügge, 2012 (Diepholz)).
For supporting such conservation activities, units for conservation management need to be delineated for the European tree frog. Most genetic analyses of the European tree frog have been conducted on a very local level and measured the genetic structure and diversity in more or less fragmented metapopulation systems (e.g. Edenhamn et al., 2000; Andersen et al., 2004; Arens et al., 2006; Angelone and Holderegger, 2009; Dubey et al., 2009). The aim of my study was to perform a large scale conservation genetic survey of the European tree frog across its distribution in Lower Saxony and adjacent areas. To allow insight on contemporary as well as historical processes I used eight microsatellite loci and mtDNA cytochrome b sequences. My specific aim was to assess genetic diversity and to define units for conservation management.
2.3 Materials and methods
2.3.1 Sample collection and preparation
Fourteen sites were sampled across the tree frog distribution in Lower Saxony and adjacent distributions in North Rhine Westphalia and Saxony Anhalt. I chose one sample site in each main occurrence of the tree frog in this region (Figure 2.1). In the occurrence near Hannover however, I sampled four sites: two in the west of Hannover (KZ, KO) and two in the east of Hannover (KH, BH) for a comparison with small scaled spatial distances. In total 237 individuals were sampled with 5 - 22 individuals per sample site (see Table 2.1). Genetic material was collected by tips of tadpole tails and by buccal swabs of adult frogs. The adults
were collected from the choruses during the breeding season in spring 2005 and 2008.
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. DNA was extracted from the buccal swabs with an Invisorb Spin Swab Kit (Invitek) following the manufacturer’s protocol and stored at -20
°C.
Figure 2.1:Current distribution of the European tree frog in Lower Saxony and adjacent areas on the basis of TK25-quadrants (grey squares) (1994-2010 in Lower Saxony (NLWKN, 2011), 1993-2006 in North Rhine Westphalia (LANUV and NRW, 2011) and 1990-2000 in Saxony Anhalt (Meyer et al., 2004)). Dashed lines denote state borders, dots denote sample sites.
Table 2.1: Overview of sample sitesa : Samples from adult frogs,t : samples from tadpoles, Ho: observed heterozygosity, He: expected heterozygosity, SD: standard deviation, FIS: inbreeding coefficient with bold values for significant difference after 1000 permutations, R: mean allelic richness over all loci, h: haplotype diversity, π: nucleotide diversity, N: number of sampled individuals, (♀): number of females included when adult frogs were sampled, NA: mean number of alleles over all loci IDSample site mean Ho ± SDmean He ± SDFISNARhπ [%] N (♀) QUQuakenbrück a 0.813 ± 0.2420.754 ± 0.137-0.0894.130.000.005 (0) WKWesterkappeln a 0.703 ± 0.2110.598 ± 0.157-0.1913.630.540.068 (0) EKEspelkamp a 0.800 ± 0.1570.764 ± 0.085-0.0505.505.360.000.0012 (0) KZKananohe Zentrum a 0.709 ± 0.1950.679 ± 0.129-0.0455.134.680.000.0020 (1) KOKananohe Ost a 0.739 ± 0.2200.692 ± 0.085-0.0724.504.46- - 11 (1) KHKolshorn a 0.778 ± 0.1390.720 ± 0.089-0.0836.505.420.530.1320 (0) BHBeinhorn a 0.731 ± 0.1360.701 ± 0.089-0.0445.634.91- - 20 (0) BABassumt 0.788 ± 0.1030.759 ± 0.090-0.0385.755.210.410.0520 RURuschwedel a 0.765 ± 0.1120.740 ± 0.070-0.0355.504.900.680.1418 (4) WGWolfsburg-Gifhorn a 0.810 ± 0.1570.800 ± 0.074-0.0137.506.520.610.0820 (1) STStrothe a/t 0.762 ± 0.1610.728 ± 0.130-0.0486.635.640.190.0221 ANAmt Neuhaus a 0.744 ± 0.1460.747 ± 0.0840.0046.505.490.660.1522 (0) SWSalzwedel t 0.631 ± 0.2280.705 ± 0.1750.1086.135.160.570.1120 PWPevestorfer Wiesen a 0.819 ± 0.1250.766 ± 0.084-0.0706.505.570.420.0520 (0)
For microsatellite analyses 5 - 22 individuals from each sample site were used. A total of eight polymorphic microsatellite loci (9, 20, 25, 60, WHA1-67, WHA1-103, WHA1-104, and WHA1-140) previously isolated by Arens et al. (2000) were amplified 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). Allele scoring was performed using the software Genetic Profiler v. 2.2.
Fragments of 901 bp of cytochrome b (cyt b) of 5 - 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 (10 mM 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).
2.3.2 Statistical analysis
2.3.2.1 Historic structure: Analysis of mtDNA
Both directions of the cytochrome b sequences were assembled using the computer software SeqMan™ II (DNASTAR, Inc., Konstanz, Germany). Multiple sequence alignments were performed in MEGA 4 (Tamura et al., 2007) 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).
Haplotype diversity (h) and nucleotide diversity (π) (Nei, 1987) were determined with ARLEQUIN v. 3.11 (Excoffier et al., 2005). 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.
2.3.2.2 Recent structure: Analysis of Microsatellites
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 and to calculate average allelic richness per population, which measures the number of alleles per locus corrected for different sample sizes. For the calculation of average allelic richness sample sites with less than ten individuals (QU and WK) were disregarded.
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. I calculated the
inbreeding coefficient FIS per sample site (Weir and Cockerham, 1984) using GENETIX v. 4.05 (Belkhir et al., 1996-2004) and tested the significance with a permutations test (1,000
permutations).
Genetic differentiation between the sample sites was calculated as pairwise FST values (Weir and Cockerham, 1984) in ARLEQUIN (Excoffier et al., 2005). However, since FST depends on marker variability and has been shown to be an imprecise estimate for genetic differentiation when applying microsatellites (Hedrick, 2005; Jost, 2008), I additionally calculated pairwise Dest (Jost, 2008), a substitute measure of genetic differentiation, using the R package DEMEtics (Gerlach et al., 2010). Significance was calculated by 10,000
bootstraps.
In order to test for isolation by distance I conducted a mantel test for correlation between pairwise genetic distances (FST and Dest) and pairwise geographic distances, implemented in IBDWS 3.23 (Jensen et al., 2005). As proposed by Rousset (1997) for populations in two-dimensional habitats, geographical distance was log-transformed and genetic distance was expressed as FST /(1 − FST) respectively Dest /(1 − Dest). Significance for r ≤ 0 was assessed via 10,000 bootstraps. Sample site QU was omitted from these analyses because of the small sample size (N = 5). The linear geographic distances among sample sites were calculated in ArcView GIS 3.3 using the Distance Matrix extension (Jenness, 2005).
Two Bayesian clustering models were conducted to infer genetic clusters. First I used STRUCTURE (Pritchard et al., 2000). 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 50,000 iterations after a burn-in period of 50,000. Because of the patchy distribution of the tree frog occurrences in this region and the large distances between most sample sites, I used the assumption of the no-admixture model and independent allele frequencies. This prior means that allele frequencies are expected in different clusters to be reasonably different from each other. Twenty runs were performed for each K. The range of possible Ks tested was from 1 to 14, 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.6.8 (Earl and vonHoldt, 2012).
Secondly, I applied GENELAND version 3.2.2 (Guillot et al., 2005a; Guillot et al., 2005b). This software uses again a Bayesian method to detect population structure but considers spatial information of the individuals. The number of genetic clusters (K) was determined by independently running the MCMC ten times, allowing K to vary from 1 to 14 to verify the consistency of the inferred K. The number of MCMC iterations was set to
100,000 per run with a thinning of 100. The uncertainty of spatial coordinates was set to 0 km and the uncorrelated frequency model was used without the assumption of null alleles.
2.3.2.3 Biogeographic zones
To identify biogeographical boundaries or zones where genetic differences between pairs of populations were largest I used the Monmonier’s algorithm as implemented in Barrier 2.2 (Manni et al., 2004). I computed the first three barriers based on cytochrome b data (p-distances) and on microsatellites (FST values). For the microsatellite data I tested the robustness of the barriers by 100 bootstrapped FST-matrices calculated via the R-package
Hierfstat (Goudet, 2005). I then considered only barriers with more than 50% support. To assess significances in this study I applied sequential Bonferroni corrections (Rice, 1989) to all multiple comparisons.
2.4 Results
2.4.1 Mitochondrial sequence analysis
I revealed 11 haplotypes of the cytochrome b fragment which differed by ten variable sites and nine parsimony informative sites (Figure 2.2; Appendix 1 and 2). Most haplotypes were unique to one sample site except haplotype Hy-1, Hy-2, and Hy-5. Hy-1 (blue) and Hy-5 (red) showed a broad distribution almost over the complete sampling area. Haplotype Hy-2 (green) was restricted to five sample sites in the north east (Figure 2.3).
P-distances were low, varying between 0 and 0.4 % (Appendix 3). The highest estimates of mtDNA diversity were found in RU (h = 0.68, π = 0.14 %) and AN (h = 0.66, π = 0.15 %). The lowest values were found in QU, EK, and KZ (all: h = 0, π = 0 %).
Figure 2.2: Haplotype network of 11 distinct haplotypes of cyt b of Hyla arborea (901 bp) in Lower Saxony and adjacent areas. Each haplotype is represented by one circle and colour. The size of the circles corresponds to the haplotype frequency. Lines between haplotypes denote mutational steps between sequences.
Figure 2.3:Distribution of cyt b haplotypes in the sample area of Lower Saxony and adjacent areas. Each haplotype is represented by one colour corresponding to the colours in the haplotype network Figure 2.2.
2.4.2 Microsatellite analysis
The eight microsatellite markers examined were polymorphic with seven to sixteen alleles per locus. The analysis using Micro-Checker uncovered signs of null alleles for the locus WHA1-67 in the sample site KO and for the locus WHA1-140 in the sample site SW. As null alleles for the two loci were found at a single sample site only, 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. The global test for HWE over all loci in each population resulted
in no significant deviation. Significant values for the inbreeding coefficient FIS were obtained for the sample sites SW (FIS = 0.108), WK (FIS = -0.191), PW (FIS = -0.070) and KH (FIS
= -0.083). No Linkage (genetic) disequilibrium was found between any pair of loci.
Since Berset-Brändli et al. (2007) found the locus WHA1-6 to be sex linked with a suppressed recombination in males, I tested its influence on the outcome of all analyses.
Excluding the sex-linked locus increased the FIS values slightly. Only sample site SW showed now significant signs for inbreeding (FIS: 0.130). Expected heterozygosity values however 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.
Genetic differentiation calculated as pairwise Dest- and pairwise FST values were significant in all cases except between the two sample sites in the West of Hannover KZ and KO (Table 2.2). In general Dest values were higher than FST values.
The Mantel test for Isolation by distance showed a significant but low correlation between genetic and geographic distances (Figure 2.4; Dest: r = 0.40, P = 0.0007; FST: r = 0.40, P = 0.0003). Indicating, that genetic differentiation is partially explained by geographic distance among sites.
Table 2.2: Pairwise Dest values (below matrix) and pairwise FST values (upper matrix) between sample sites; bold = significant difference after sequential Bonferroni correction. QUWKEKKZKOKHBHBARUWGSTANSWPW QU00.1140.0530.1210.1250.1170.1100.0460.1120.0540.0870.0910.0910.067 WK0.19800.1080.0870.1210.1270.1630.1040.1800.1130.0890.1340.1540.153 EK0.1360.28500.0510.0430.0750.0850.0650.0850.0680.0930.0710.1000.074 KZ0.3510.2230.1470-0.0010.0870.1050.0750.1250.0880.0880.0970.0960.096 KO0.3870.3020.145-0.00400.0850.1070.0730.1180.0880.0870.0910.0910.096 KH0.4490.3480.2550.2740.27000.0360.0600.0630.0720.0580.0940.1080.086 BH0.4190.4260.2750.3130.3250.09600.0810.0920.0800.0750.1170.1180.092 BA0.2190.2800.2570.2480.2520.2090.26600.0680.0370.0710.0790.0780.070 RU0.4260.4640.2920.3560.3610.1940.2770.22900.0600.0670.0730.0990.068 WG0.1880.3460.3070.3310.3340.2700.2840.1500.22600.0500.0730.0820.063 ST0.3270.2550.3460.2870.3010.1960.2450.2560.1900.18300.0620.1020.068 AN0.2960.3720.2570.3030.3120.3120.3710.2820.2500.2930.23300.1070.084 SW0.2960.3970.3350.2590.2480.3760.3860.2660.3230.2870.3600.38700.079 PW0.2730.4130.2990.2570.2920.3200.3100.2760.2580.2590.2330.3050.2670
Figure 2.4:Isolation by distance plots. (a) Dest /(1 - Dest) versus log geographic distance and (b) FST /(1 - FST) versus Log geographic distance. Lines are the RMA (reduced major axis) regression.
For the STRUCTURE analysis both approaches to determine the correct number of genetic clusters, the ΔK statistics by Evanno (2005) and the average log likelihood- values Pr(X|K), peaked clearly at K = 7 (Figure 2.5; Figure 2.6). One main cluster in the West of Hannover consisted of the sample sites WK, EK, KZ, KO (red). The sample sites KH and BH in the East of Hannover formed an extra cluster (blue) separated from the sites in the West. Sample site WG was admixed with large parts of the cluster build by BA (yellow) and RU and ST (green).
The sample sites in the East – SW, PW and AN were separated in three further clusters (orange, grey and light-blue) whereas the latter formed one cluster together with QU in the West of Lower Saxony.
Figure 2.5:Mean values of estimated Ln probability of data (LnPD) for each K (a) and delta K (b)
Figure 2.6: STRUCTURE bar plot for K = 7; QU, WK, EK etc. = sample sites, separated by fine black lines. Each individual is represented by a single vertical line broken into K-coloured segments, with lengths proportional to each of the K-inferred clusters.
In consistency with the STRUCTURE result, in the GENELAND analysis the highest average log posterior probability was found for seven genetic clusters. Cluster 1: WK, EK, KZ, KO;
cluster 2: KH, BH; cluster 3: QU, BA, WG; cluster 4: RU, ST; cluster 5: SW; cluster 6: AN;
and cluster 7: PW. This is the same clustering as found by the Structure analysis except for the assignment of sample site QU. (Figure 2.7)
Figure 2.7: Map of estimated cluster membership for K = 7 by GENELAND analysis. Each cluster is shaded in a single colour.
2.4.3 Biogeographic zones
The most significant genetic discontinuities among sampled locations were estimated on the basis of microsatellite FST-values and cyt b -p-distances. Obtained barriers are marked with lines in figure 2.8. Two barriers were supported by both markers (microsatellites and cyt b).
One barrier separated KZ from KH (KO and BH were not considered in sequence analysis), and a second barrier was found between WK and EK.
Figure 2.8: Most important barriers to gene flow from the BARRIER analysis (Manni et al., 2004). Red Lines indicate most significant barriers to gene flow (> 50 % bootstrap support) estimated on microsatellite FST-values.
Significance of barriers are given as percent bootstrap. Barriers estimated on the basis of cyt b p-distances are marked with grey lines ranked I-III, in order of decreasing magnitude. In green the Delaunay triangulation, in blue the Voroni tessellation between sample sites, used to calculate borders to neighbouring sample sites.
2.5 Discussion
I analysed the broad scaled genetic structure and variation of the Lower Saxonian tree frog occurrences to aid conservation management implementation. Cyt b sequences showed low variation but distinct geographic-genetic pattern was revealed. Using microsatellite analysis I found seven distinct genetic clusters. Genetic diversity was high in most sample sites.
For supporting measures of effective conservation management, the identification of conservation units or management units (MUs) is critical (Palsbøll et al., 2007).
“Management Units are defined as demographically independent breeding units and are identified as populations having distinctive allele frequencies” (Moritz, 1994b). Their recognition is fundamental to proper short-term management (Moritz, 1994a). The genetic clusters revealed by STRUCTURE and GENELAND can more or less been used to delineate such
management units as conducted e.g. for Koalas in the area of Sidney (Lee et al., 2010).
Palsbøll et al. (2007) propose that “MU status should only be assigned when the observed estimate of genetic divergence is significantly greater than a predefined threshold value.”
Unfortunately, corresponding studies are missing especially for amphibians defining such a
Unfortunately, corresponding studies are missing especially for amphibians defining such a