Genetic analyses of brown hare (Lepus europaeus) support limited migration and translocation of Greek populations.

Numerous studies have shown that the phylogeography of many species, including European brown hare, has been affected by the climatic oscillations of the Pleistocene. During this period the Balkans acted as a major refugium offering habitable conditions for many species. However, few studies have focused on the specific role of the Greek peninsula in the phylogeographic history of species in this southernmost margin of Balkans. We, therefore analyzed a 528 bp fragment of the D-loop region of mitochondrial DNA in 154 wild brown hare individuals from unsampled areas from both mainland and island Greece and compared it to 310 available brown hare sequences (including 110 Greek samples). Newly identified haplotypes show characteristic distribution in specific Greek areas reinforcing the theory that Greece can be considered as a subrefuge within Balkans for a number of species, with several "refugia within refugia" spots, holding significant genetic diversity. No haplotypes from wild Greek individuals clustered with the Central and Western Europe group revealing a minimal contribution of this area to the colonization of central Europe. One hundred and ten reared brown hares were also analyzed to elucidate the impact of introductions on local populations. Most of these samples presented close genetic affinity with haplotypes from Central and Western Europe indicating that farms in Greece use breeders from those areas. Therefore, despite human translocation of individuals, the genetic structure of brown hare has mostly been influenced by paleoclimatic conditions and minimally by human actions.

Introduction

Pleistocene glaciation cycles have considerably influenced the current geographic distribution of many European temperate species. During the Last Glacial Maximum (LGM, 26–19 ka, [1]), large parts of Europe presented inhospitable habitats and diverse species were limited to warmer regions of Europe [2, 3]; the southern peninsulas of Italy, Iberia and Balkans (extending into Asia Minor) provided more habitable conditions and acted as major refugia for many species [4-6]. The populations that survived in these isolated regions underwent limited gene exchanges resulting in independent evolution and genetic divergence. According to this theory, these separate refugia contributed significantly to the post-glacial recolonization of the rest of Europe [2, 7, 8]. More recent evidence proves that many species also survived in unexpected latitudes in northern or cryptic refugia [3, 9–11]. Such refugia have been identified in regions of Western Europe, Central Europe and the Carpathians [8, 12–14]. “Refugia within refugia” have been identified even within the three Mediterranean peninsulas for animal as well as plant species [15-18].

Pleistocene climate changes are believed to have, also, profoundly affected the phylogeographic structure and the evolutionary history of European brown hare (Lepus europaeus). This mammal presents considerable plasticity evidenced by its wide distribution in a range of environments, while it has also been effectively introduced, even in tropical habitats [19]. Its mobility can differ between the two genders (females present a more philopatric behavior, [20]) and between different periods of the year (higher levels of mobility in breeding seasons, [21]). According to previous studies the current genetic structure of European populations is the result of its postglacial spread from the Balkans and Asia Minor [22-24]. Both Kasapidis et al. [22] and Stamatis et al. [24] concluded that brown hares from Greece, Asia Minor and Central Europe were separated in two distinct clades, which may correspond to the post-glacial refugia of Balkans and Anatolia. Fickel et al. [23] confirmed these separate clades and identified an additional Italian one.

However, data on the evolutionary history of brown hare in Greece is missing, especially on populations from Central Greece and the Greek islands, which, some of them, have been shown to be connected to mainland Greece due to sea oscillations during glaciations cycles [25]. We analyzed a fragment of the mitochondrial control region sequence of brown hares both of mainland and island Greece and reared brown hares as well, and we combined our data with results from previous studies aiming to three main objectives: (a) highlight the genetic diversity of brown hare in Greece, (b) search for signs of unique/genetically differentiated lineages in agreement with the ‘refugia within refugia’ hypothesis by applying a more extensive sampling of these regions (c) investigate the impact of introduced genomes by reared individuals into the local populations as it known that restocking programs have taken place in many European countries in order to mitigate brown hare population reductions, without always taking into consideration the effect on the local genetic pool [26-28],

Materials and methods

Ethics statement

No animals were sacrificed only for the purposes of this study. The wild brown hare samples in this study represent material collected opportunistically (after hunting procedure) from animals’ hunter-harvested by members of the Greek Hunting Federation of Sterea Hellas during the hunting seasons, according to the prerequisites of the Greek Legislation (Ministerial Decision 103305/6093/14-11-2007, FEK 1626/Β/13.8.2008, FEK 1611/Β/2009, FEK 1183/6-8-2010, FEK 1763/4-8-2011). The reared samples included in this study were collected by private veterinarians. No domesticated individual was euthanized during the study and efforts were taken to mitigate animal suffering. There is not an ethics committee in the Aristotle University of Thessaloniki. But, all experimental procedures follow the code of conduct of the Research Committee of the Aristotle University of Thessaloniki. The study was approved through the acceptance of the program by the Research Committee of the Aristotle University of Thessaloniki (Program number: 84430).

Sample collection and DNA extraction

A total of 154 wild brown hares were collected from various regions of both mainland and island Greece with the assistance of the 4th Hunting Federation of Sterea Hellas. Muscle tissues were obtained during the hunting seasons of 2006 to 2009. Particularly, we conducted a more comprehensive sampling in Central Greece and Aegean islands (including areas which had not been previously sampled or had been poorly sampled) (S1 Table). Additionally, 110 reared individuals were analyzed from four Greek breeding stations (two breeding station from Evia, one from Agrinio and one from Macedonia) (a sloughed piece from animals’ ear was used). Whole genomic DNA was extracted using the protocol of [29] which is based on the chemical compound CTAB.

PCR amplification and DNA sequencing

A fragment of 528 bp of D-loop (control region, CR) of mitochondrial (mt) DNA (including the entire CR-1) was amplified using primers Le.H-Dloop and Le.L-Dloop [22]. The 5’ ends of the two primers correspond to positions 15423 and 15951 of the complete mtDNA, respectively. The total volume of polymerase chain reaction was 25μl in which 100ng of genomic DNA was amplified, using 0.05 units of Qiagen Taq polymerase, 2mM dNTPs, 5 pmol of each primer, 2.5 mM MgCl2 and 2.5 μl of 10 X Reaction Buffer. Thermal cycling amplification conditions were as follows: initial denaturation at 94°C for 3 min, followed by 33 cycles of strand denaturation at 94°C for 1 min, annealing at 50°C for 45 s and primer extension at 72°C for 40 s and a final 3min elongation time at 72°C. The PCR products were purified using the Nucleospin Extra kit (Macherey-Nagel, Duren, Germany) and sequenced by Sanger method from Macrogen Inc. (Seoul, Korea) using an ABI 3730XL DNA Analyzer.

Sequence analyses

A total of 264 sequences (dataset 1) (Genbank Accession nos: MH842011-MH842084) were aligned using ClustalW [30] with final adjustments by eye. A second dataset was created aiming at a pan-European analysis. Thus, data from the present analysis were combined with 310 CR-1 sequences of L. europaeus retrieved from GenBank: 72 sequences described in Kasapidis et al. [22], 69 sequences described in Stamatis et al. [24], 109 sequences described in Fickel et al. (direct submission), 42 sequences described in Fickel et al. [23] and 18 sequences described in Antoniou et al. [31]. All above sequences were long enough to retain informative sites for proper downstream analyses. The overlapping region with the new dataset was 335bp.

The best-fit nucleotide substitution model for this dataset was determined using jModelTest 2.1.10 (under the Bayesian Information Criterion, BIC) [32, 33]. The GTR + I + G model [34, 35], was used for Bayesian phylogeny analyses, carried out with Beast 1.8.4 [36]. The Bayesian tree was constructed using a strict clock model and a coalescent tree prior. The analysis was run for 108 Markov Chain Monte Carlo (MCMC) generations, sampled every 104 generations. Convergence of chains was visualized using Tracer 1.6 [37] discarding the first 20% of trees as burn-in. ESS values for all parameters were > 248, larger than the threshold value of 200 identified by Tracer v. 1.6. The trees produced by Beast were then summarized in TreeAnnotator 1.8.4 and visualized in FigTree 1.4.3 [38]. A D-loop sequence of Lepus timidus (Accession number: EF515861), was used as outgroup.

Additionally, a median-joining network [39] was constructed using the software Network 5.0.0 and the frequencies of the sequences (Fluxus Technology) assuming equal weights for all mutations and setting the genetic distance parameter e to zero in order to restrict the choice of feasible links in the final network.

Genetic variability indices of brown hare populations, i.e. the number of distinct haplotypes (h), haplotype (Hd) and nucleotide diversity (π) values, were estimated using DnaSP 5 [40]. To estimate the haplotype richness by additional extensive sampling in Greece, a rarefaction curve was calculated in RarefactWin [41] and plotted in Microsoft Excel (Office 365).

Results

A total of 74 different haplotypes were identified among the 264 newly analyzed samples, 82.4% of these have never been detected before. Fifty seven haplotypes were found in wild brown hares and 20 in domestic individuals; only three haplotypes were shared between wild and domestic samples. In this dataset 1, the population of Central-South Greece possesses the highest haplotypic diversity (0.938), whereas North Greece possesses the highest nucleotide diversity values (0.038) (Table 1). Most haplotypes from Aegean islands (9 out of 12) were unique, characterizing each studied island. Additionally, North Greece populations revealed high percentage of (population-specific) haplotypes per number of individuals. As regards reared populations, haplotype and nucleotide diversity values were generally lower compared to the wild ones (Table 1) with breeding station 2 exhibiting remarkably low diversity values (Hd = 0.275, π = 0.004).

Table 1

Genetic diversity values for wild brown hares from five Greek areas and Central/West European countries and for domestic brown hares from four Greek breeding stations.

Dataset 1 includes new samples from this study, whereas dataset 2 includes additional sequences from GenBank (for references see MM section).

Region	Dataset	N	H	Hp	V (%)	W (%)	Hd (SD)	π (SD)
North Greece	1	16	11	9	68.7	56.2	0.900 (0.062)	0.038 (0.003)
	2	92	66	64	71.7	69.6	0.982 (0.007)	0.034 (0.002)
Central–South Greece	1	89	37	31	41.6	34.8	0.938 (0.016)	0.018 (0.002)
	2	117	60	56	51.3	47.9	0.963 (0.010)	0.022 (0.001)
Central Aegean Islands	1	45	11	8	24.4	17.8	0.855 (0.031)	0.012 (0.002)
	2	61	20	17	32.8	27.9	0.914 (0.019)	0.015 (0.001)
Eastern Aegean Islands	1	1	1	1	100.0	100.0	-	-
	2	12	6	6	50.0	50.0	0.864 (0.064)	0.023 (0.003)
Ionian Islands	1	3	1	1	33.3	33.3	-	-
	2	5	1	1	20.0	20.0	-	-
Total of wild Greek samples	1	154	57	54	37.0	35.1	0.965 (0.006)	0.025 (0.002)
	2	287	149	149	51.9	51.9	0.982 (0.003)	0.030 (0.001)
Central/West Europe	2	841	130	130	15.5	15.5	0.893 (0.009)	0.009 (0.001)
Breeding Station 1	1	19	5	2	26.3	10.5	0.801 (0.046)	0.007 (0.002)
Breeding Station 2	1	60	6	4	10.0	6.7	0.275 (0.074)	0.004 (0.002)
Breeding Station 3	1	20	9	4	45.0	20.0	0.868 (0.050)	0.019 (0.001)
Breeding Station 4	1	11	4	4	36.4	36.4	0.745 (0.098)	0.007 (0.001)
Total of reared samples	1	110	20	17	18.2	15.5	0.704 (0.043)	0.013 (0.001)

n, sample size; h, number of different haplotypes; Hp, population/region specific haplotypes; V, the percentage of the total number of haplotypes/number of individuals; W, the percentage of the population–specific haplotypes/number of individuals; Hd, haplotype diversity; π, nucleotide diversity (SD, Standard deviation)

The combination of our haplotypes and the retrieved sequences (384 sequences in total, corresponding to 1128 individuals) resulted in a dataset of 326 unique haplotypes (dataset 2). In this dataset, Greece possessed higher haplotypic (0.982) and remarkably higher nucleotide (0.030) diversity values compared to Central/West Europe. It is remarkable that Greece possesses slightly higher number of haplotypes than Central/West Europe (149 vs 130), despite the fact that its sample size is almost one third of Central/West Europe (287 vs 841 samples). In addition to this, according to the rarefaction curve the plateau has not been reached for the discovery rates of new haplotypes despite the additional extensive sampling in Greece, confirming the great diversity of Greek populations (S1 Fig). Furthermore, Greece possessed remarkably higher percentage of (population–specific) haplotypes per number of individuals compared to the rest of Europe.

The Bayesian phylogenetic analysis (Fig 1) supported haplotype clustering in two distinct clades, a result consistent with previous studies [22, 24]. In the first clade, conventionally named “Anatolian” (A), haplotypes from Turkey, Cyprus, Israel, Bulgaria, Greece as well as an Italian haplotype were grouped. Only 4 haplotypes (two of them newly characterized) from the present study clustered in this clade (samples from Chios and Evros).

Fig 1

Phylogenetic tree (A) and distribution (B) of brown hare haplotypes in the Balkans. Subgroup A consists the “Anatolian” clade, while subgroups CWE and E1-E5 consist the “European” one. Each coloured circle in b represents the existence of individuals belonging to the coloured branches of the tree in a. The size of circles in b reveals the number of individuals studied from each area. Posterior probabilities >0.7 are shown (Map was modified from https://commons.wikimedia.org/wiki/File:Blanc_neutral_map_of_Greece.png, CC BY 2.5).

The second clade, the “European” one involved subgroups CWE and E1-E5 (Fig 1). This clustering is supported by high values of posterior probabilities (posterior probabilities >0.7 are shown, Fig 2). Subgroups E1-E5 were found mainly in Greece and Bulgaria. Most of the haplogroups presented specific geographic distribution within Greece. The E1 subgroup is mainly found in Evia island and Central Greece, E2 is spread over mainland Greece, the E3 mainly in North Greece, E4 mainly in Aegean islands (Cyclades and Crete), Central and South Greece. E5 is more widespread and observed in Northeastern Greece, Slovakia and Italy. The CWE subgroup includes the vast majority (122 out of 130) of haplotypes from countries of Central and Western Europe. Only 8 haplotypes from these areas were clustered in the E2 (one haplotype) and E5 (7 haplotypes) subgroups. No wild Greek haplotypes belong to subgroup CWE. The majority of reared individuals (97 out of 110 individuals, 13 out of 20 haplotypes) clustered in the CWE subgroup, while a low number of domestic haplotypes was included in groups E1, E2, E4 and E5 (7 haplotypes found only in 13 individuals).

Fig 2

Median-Joining network of 335bp mtDNA brown hare haplotypes (326 haplotypes) in the full dataset.

The proportional size of nodes indicates the frequency of haplotypes. The coloured line of circles corresponds to the coloured subgroups in the phylogenetic tree (Fig 1). Small red dots represent inferred haplotypes.

In addition to the Bayesian analysis, a median-joining network was constructed to further elucidate the phylogenetic relationships among haplotypes (Fig 2). The grouping of haplotypes in the MJ network is in concordance with the clustering in Bayesian analysis. The network verified the existence of two distinct clades (Anatolian and European), which were separated by several mutational steps.

Discussion

The combination of our data with previous studies has resulted in a broader and clearer depiction of brown hare genetic diversity distribution and sheds more light on the evolutionary factors that have shaped it.

Importance of rear-edge Greek brown hare populations

Our data reveal that Greek populations present remarkably high genetic diversity. Although other studies have also illustrated Greece’s role as a biodiversity hotspot [2, 5, 17, 42, 43] it is the first time that such high diversity values are reported. The identification of numerous newly characterized haplotypes is important for revealing the evolutionary history of brown hare in Greece and Europe. These higher haplotype and nucleotide diversity values of Greek population as well as the high percentage of (population–specific) haplotypes per number of individuals, underline the existence of long-term favorable environmental conditions, under which these populations remained relatively stable and genetically variable [44]. The existence of a mosaic steppe in Greece [45, 46], has provided the appropriate environment for brown hares to survive and evolve.

The high diversity values in the Greek marginal populations combined with the low diversity of central European ones do not support the “centre-periphery hypothesis” which includes the idea that core populations present higher genetic diversity than the edge ones [47, 48, 49]. This pattern has also been observed in other species in Europe like the obligate myrmecophilus butterfly (Maculinea arion, [50]), which presented higher genetic diversity in one of the main European refugia, Italy. Such rear edge populations, manage to survive under suitable for species persistence conditions [51], maintaining a large fraction of their pre-glacial diversity and thus leading to the creation of independent refugia within the known refugia of southern European peninsulas [52]. This scenario has been demonstrated for a number of Greek species including wild boar [17], horned viper [43], alpine newt [42] and red deer [53] that have been reported to harbor genetically differentiated lineages from the rest of the Balkans.

The fact that high genetic diversity values of Greek populations were maintained for long time periods is also very well presented by the complex structure of the Bayesian phylogenetic analysis and the median-joining network, contradicting undeniably the “centre-periphery hypothesis”. The star-like phylogeny in CWE subgroup combined with the absence of geographic structure support the idea of rapid expansion and colonization of brown hare in Central and West Europe. On the contrary, the Greek haplotypes (mainly included in subgroups E1-E5) do not present very close phylogenetic relationships and there is a clear geographic patterning.

Although the concepts of geographic centrality, ecological marginality and centre of origin have been repeatedly presented in the literature, very few studies have focused on all these factors at the same time [44]. According to the centre-periphery hypothesis [54], species are more abundant in the centre of their geographic ranges, where climatic conditions are more optimal for their survival. As distance from the centre increases the environmental conditions become unsuitable. Subsequently, marginal populations may be more vulnerable to any changes and consequently more prone to extinction. Our results, however, support the view of Pironon et al. [44] who argue that it may, actually, be more accurate if we consider refugial populations (such as those from the southern Balkans, which are at a geographic rear-edge) as the central ones and populations from recently colonized areas (geographically central) as peripheral, and that therefore, the centre and the margin of populations cannot always be geographically delimited.

Microgeographic structure of brown hares in Greece

The identification of haplotypes in newly sampled areas also gives us a more comprehensive view of the phylogeography of the species within Greece. As shown both by the Bayesian and network analyses different lineages (E1, E2, E3, E4) exist, with characteristic geographic structure. Similar results were obtained by Kasapidis et al. [22], who presented different clusters with specific geographic distribution. However, the present study reveals a much more extended geographic structure as we added a significant number of new samples and the clustering is strongly supported by high posterior probability values. These lineages may indicate the existence of different small refugia in Greece during the Last Glacial Maximum, supporting the “refugia within refugia” model. This model has been proven to fit many other species as well, for all three main refugia (Balkans, Iberia, Italy) [17, 43, 55].

Apart from historical interplay effects on the genetic structure of the species some other present day factors have to be considered in order to be more accurate in explaining current genetic data. Specifically, the present day isolation of different Greek populations cannot be only attributed to the geographical boundaries as there are no natural barriers between all sampling areas. According to [20], it can be explained by philopatric females that do not disperse in large distances. The maternal inheritance of mtDNA, used in the present study could overestimate the above differentiation. However, results of Antoniou et al. [31] are consistent in mitochondrial as well as nuclear markers, indicating that the philopatry of females does not affect the genetic structure of the species. Nevertheless, the study of Antoniou et al. [31] focused on a restricted area in Northeastern Greece and therefore a more thorough overview of demographic genetic variability in an extended region and the use of both mitochondrial and nuclear markers are essential, for safe conclusions.

Despite the high genetic variability of Greek populations, Aegean islands presented low haplotypic diversity values (Table 1), probably due to low population densities and thus reduced effective population sizes. It is noteworthy to mention that every different island studied (Andros, Paros, Tinos and Kythira) is represented by a unique haplotype, which is characteristic for this specific island, regardless of the number of samples analyzed. The long term geographic isolation of these populations, and related absence of gene flow could explain the presence of these unique haplotypes. One exception was noticed in the case of the island Naxos, in which two of the three haplotypes were shared with individuals from Evia. A contact of these areas in previous glacial periods is possible [56]. The low density of island populations and their low genetic variability render them more vulnerable to the introduction of foreign genomes in these areas, which can lead to deterioration of their genetic pool in short periods [57]. For this reason, it is important to protect them in the context of the protection of biodiversity of the species.

Considering the clustering of Chios with the Anatolian clade, it could be attributed to: i) human induced translocations of individuals from Turkey on this island and/or ii) vicariance when Chios separated from Turkey (about 9,000 years ago, [58]). However, any conclusion cannot be reached safely ​​without additional sampling and further analysis including archeological samples. This also applies for previously sampled populations from the islands of Lesvos and Samos, which also cluster with the Anatolian clade. Rhodes was not connected to the Asia Minor mainland since the end of the Pliocene or the beginning of the Pleistocene [59]. The presence of brown hare there is most probably due to human translocation.

The fact that many unique haplotypes have been found only in Greece, absent from Central Europe, supports the hypothesis that Greek populations were not involved in the northward expansion of the species. These data reinforce the concept that colonization of Europe took place starting from Central and Northern Balkan refugia [5], while the contribution of southern distinct populations was minimal. This scenario has been also proposed for species such as the grasshopper and the hedgehog and one possible explanation has been given from [5], indicating that northern Balkan populations of brown hare prevented the expansion of Greek ones from the southernmost distribution of the species.

Our data indicate that the north part of Greece was the region with the highest mtDNA diversity values (Table 1). This could be attributed to the coexistence of the two highly differentiated clades in the area of Thrace (European and Anatolian clade). According to previous studies [5], a well-known hybrid zone exists in North-Eastern Greece (Thrace) and South-Eastern Bulgaria/North-West Turkey [22, 24, 31] where the eastern (Anatolian) populations have managed to introgress. This also holds for other taxa (grasshopper, [56]) but it is not always consistent; for species including wild boar [17] and crested newt [60] it is clear that there is no such geographical overlap, whereas for the European ground squirrel, the overlap zone is limited [61]. It is clear that there is no common rule about the existence of an overlapping hybrid zone for all studied species. More research is needed in a variety of species in order to elucidate the mechanisms and the biological history of the species that can produce such differences, despite the shared habitat.

Impact of restocking programs

During the last decades numerous restocking programs have taken place introducing foreign individuals of brown hare into several European countries [26, 28, 62]. In Greece, individuals have originated from breeding stations in Italy, Yugoslavia and Bulgaria [63]. These restocking programs were uncontrolled, in many cases, with unpredictable consequences on the historical distribution and genetic integrity of indigenous hare species [26-28]. The introduction of foreign genomes into local populations may influence its fitness [27].

The results of the present study indicate that the majority of reared individuals cluster with haplotypes from Central and Western Europe. Since release of such individuals in the environment is a common practice, this constitutes a potential risk to the domestic gene pool of the species. This is aggravated by the fact that farmed individuals can survive in the wild at least for a reproductive period, so it is possible that they transfer their genome [63]. However, it is also true that no Greek wild samples were found to belong to the CWE group, even though most reared individuals carry such haplotypes. This supports that no gene flow has yet been detected between reared and wild Greek populations. On the other hand, some haplotypes from domestic brown hares clustered with Greek wild haplotypes, implying that these individuals were collected from the wild and used as breeders. In any case, the use of endemic individuals as brood stock is a practice that should be followed in order to protect region specific diversity [28].

Conclusion

Although during LGM, the Balkan peninsula acted as a major refugium for the colonization of the northern parts of Europe by brown hare [3, 5], the contribution of Greece seems to be minimal. Greek populations have retained a large fraction of their pre-glacial diversity and though restricted to a small refugial area, have managed to survive. Our analysis shows that a “refugia within refugia” model fits better brown hare evolution. These “peripheral” populations hold important abundant genetic variability relative to central populations and constitute significant units that deserve conservation priority [47]. For this reason, it is essential that any restocking program should be accompanied by genetic control of the individuals that will be released in the wild, in order to protect this genetic stock. This is indeed requested by the Greek state (Ministry Decision 98161/4136/29-9-2008 published in FEK 637/B/6-4-2009) for all stocked game animals.

Names of sampling localities for each region and corresponding geographic coordinates, number of samples and number of haplotypes.

(DOCX)

Click here for additional data file.

Rarefaction curve for haplotype diversity of Greek brown hare individuals.

(DOCX)

Click here for additional data file.