A European Concern? Genetic Structure and Expansion of Golden Jackals (Canis aureus) in Europe and the Caucasus.

In the first continent-wide study of the golden jackal (Canis aureus), we characterised its population genetic structure and attempted to identify the origin of European populations. This provided a unique insight into genetic characteristics of a native carnivore population with rapid large-scale expansion. We analysed 15 microsatellite markers and a 406 base-pair fragment of the mitochondrial control region. Bayesian-based and principal components methods were applied to evaluate whether the geographical grouping of samples corresponded with genetic groups. Our analysis revealed low levels of genetic diversity, reflecting the unique history of the golden jackal among Europe's native carnivores. The results suggest ongoing gene flow between south-eastern Europe and the Caucasus, with both contributing to the Baltic population, which appeared only recently. The population from the Peloponnese Peninsula in southern Greece forms a common genetic cluster with samples from south-eastern Europe (ΔK approach in STRUCTURE, Principal Components Analysis [PCA]), although the results based on BAPS and the estimated likelihood in STRUCTURE indicate that Peloponnesian jackals may represent a distinct population. Moreover, analyses of population structure also suggest either genetic distinctiveness of the island population from Samos near the coast of Asia Minor (BAPS, most STRUCTURE, PCA), or possibly its connection with the Caucasus population (one analysis in STRUCTURE). We speculate from our results that ancient Mediterranean jackal populations have persisted to the present day, and have merged with jackals colonising from Asia. These data also suggest that new populations of the golden jackal may be founded by long-distance dispersal, and thus should not be treated as an invasive alien species, i.e. an organism that is "non-native to an ecosystem, and which may cause economic or environmental harm or adversely affect human health". These insights into the genetic structure and ancestry of Baltic jackals have important implications for management and conservation of jackals in Europe. The golden jackal is listed as an Annex V species in the EU Habitats Directive and as such, considering also the results presented here, should be legally protected in all EU member states.

Introduction

An implementation of molecular techniques to study population genetics has broadened our knowledge about several aspects of wildlife biology and ecology, including breeding characteristics [1, 2], population connectivity, and dispersal [3, 4]. Simultaneously, it enabled us to assess effects of historical processes [5-7], habitat fragmentation and isolation on distribution of genetic diversity (e.g. [8-10]) and to reconstruct routes of recent colonisations, range expansions and biological invasions [11-14]. As a result, information provided from molecular markers is frequently used in wildlife management and conservation of endangered species [15-19].

Changes in the geographical range are recognized as natural processes and have occurred in the history of most species [20-23]. Theoretical aspects of genetic after-effects of range shifts have been thoroughly analysed (e.g. [24]). It was shown that range expansions may lead to changes in population genetic structure and diversity. Initially, genetic structure should be clearly emphasized and genetic diversity in subdivided population will likely be reduced in comparison with the main distribution range and/or source population due to repeated bottlenecks. However, over time as new areas are occupied, connectivity among territories may be established and spatial population structure might decrease due to balanced gene flow among populations, causing homogenization and increased genetic diversity within populations [24-26]. Surprisingly, the genetic consequences of natural, contemporary range expansions have begun to be investigated only very recently [27-31] and results so far are equivocal and not always concordant with theoretical expectations.

Several carnivore species are currently expanding their distributions, especially in Europe [32]. It has been observed that such populations are characterized by particular genetic structure and processes, at least on the scale of individual countries. The study of the recently expanding (most probably from Russia) brown bear (Ursus arctos) population in Finland revealed disappearance of initial structuring and homogenization, as well as gradual increase of genetic diversity [33], as expected on the basis of theoretical models of range expansion [24, 25, 34]. Moreover, Hagen et al. [33] have shown increasing admixture between two genetic clusters occurring in Finland [35] as the range expansion proceeded. In contrast, the Finnish grey wolf (Canis lupus) population, also expanding since the 1990’s after almost complete eradication in the 19th century, exhibited decreased genetic diversity during the initial phase of expansion, despite clearly lower estimated population size [36]. The authors attributed this result to a low degree of connectivity with adjacent Russian wolf population.

The golden jackal (Canis aureus) is one of the most widely-distributed canid species, found in many areas of Europe and southern Asia [37, 38]. The ongoing expansion of the species in Europe has caused concerns in regard to possible negative effects its presence could exert, for example through excessive predation of other wildlife species or livestock, and the transmission of pathogens. In addition, there are several uncertainties regarding jackal management and policies, often in association with the unknown origins of jackal populations [38].

Population genetics of this species has been so far poorly characterised, especially when compared to Europe’s large carnivores, such as the grey wolf (e.g. [39-42]), the European lynx (Lynx lynx) [43-47] or the brown bear (e.g. [35, 48–50]). The first study focused on jackals in Serbia [51] suggested a low level of genetic diversity and weakly pronounced genetic structure in this recently-spreading population (see also [52]). Low genetic differentiation was also found in populations from Bulgaria, Croatia, and Italy [52]. A significant but weakly-pronounced genetic structure was only observed in the population of jackals from Dalmatia (Adriatic coast of Croatia). Fabbri et al. [52] also discovered that the jackals in Italy have an admixed origin from the Dalmatian and mainland populations. The genetic data in these cases were suggestive of a colonization process in golden jackals that is predominantly of a ‘stepping-stone’ nature, with short-distance dispersal and intermediate admixture. This contrasts with the long-distance dispersal observed in other canids, such as grey wolves [53, 54].

Genetic relationships of the European golden jackals with jackals from the Asiatic part of the species’ range, were not yet determined. Moreover, none of the studies so far analysed genetic structure of the population on the larger scale (i.e., the continental level). Consequently, the understanding of historic development of jackal populations in Europe is lacking. One of the hypotheses suggested that the European population goes back to the introduction of jackals from northern Africa in the 15th century [55]. This was later rejected on the basis of morphology [56, 57], but the origin of most of the European population remains unknown. Archaeologic data indicate that jackals were already present along the Mediterranean coast in Croatia and Greece ca. 7,000–6,500 yBP [58, 59]. Jackals remained absent from most of Europe until the 19th century, when the species started to expand slowly, followed by a rapid expansion at the end of the 20th century, which continues today [38, 60]. However, it is unclear whether any of the present European populations originate from this ancient Mediteranean population or if they are decendants of the later Asian colonization, e.g. from the Middle East or the Caucasus. Secondly, if there was a recent colonization from the east, it is unknown whether original small Mediterranean populations survived and merged with the wave of recent expansion. It is also unknown whether low genetic diversity and lack of distinct genetic structure in part of the European golden jackal population [51, 52] is an after-effect of fragmentation and population decline in the first half of 20th century, or rather resulted from recent expansion, interlinked with the founder effect pertaining to a recently established population. Hence, samples from potentially long-lasting, stable populations, such as southern Greece, should be analyzed. Although it was suggested that Italy was colonised from the Dalmatian coast and the mainland [52], the source of other expansions in Europe have not yet been identified. The lack of proper knowledge about the history of golden jackals in Europe can significantly affect management decisions and thus influence the conservation of the species. For example, the Estonian, Latvian and Lithuanian governments, despite the lack of reliable data, consider the golden jackal to be an alien species introduced to the Baltics by people, and based on this, these governments recently allowed unlimited lethal removal with the goal of eradicating the species [38].

The aim of the present work is to characterise for the first time the population genetic structure of European golden jackals on the continental scale, with the incorporation of samples from hitherto unstudied regions. Therefore, we included samples from the Peloponnesus Peninsula (southern Greece), which could possibly originate from the Neolithic population [59]; the insular population on the island of Samos located 1.7 km from the coast of Asia Minor, which represents the first investigation of an island population of the species; and the population from the Caucasus, a region known as a ‘hotspot’ for biodiversity [61]. An attempt is also made to identify the origin of the recently-established population in the Baltic States, and hence to resolve its controversial status and aid management decisions.

Material and Methods

Samples

Tissue samples used in this study were obtained from 97 individuals originating from five geographical regions (Fig 1, Table A in S1 File), i.e. i) south-eastern Europe (SEE)—comprising samples from Romania (country code ROU; n = 5), Croatia (HRV; n = 2), Slovenia (SVN; n = 2), Ukraine (UKR; n = 12), Serbia (SRB; n = 25), Hungary (HUN; n = 10), and northern Greece (GRC; Chalkidiki Peninsula, n = 1); ii) the Caucasus (CAU)—comprising samples from Mountainous Karabakh (NKR; n = 6), Armenia (ARM; n = 3), and Georgia (GEO; n = 5); iii) the Baltic States (BAL)—comprising samples from Lithuania (LTU; n = 1) and Estonia (EST; n = 4); iv) southern Greece (GRE-P) —comprising samples from the Peloponnese (GRC; n = 11); and v) the island of Samos (GRE-S, n = 10).

Fig 1

Distribution of sampling sites. Shaded areas represent areas with permanent presence of jackals (based on [38] and [37]).

Molecular genetics protocols

Genomic DNA was extracted using NucleoSpin Tissue kit (MACHEREY-NAGEL) with standard protocol. We amplified 15 microsatellite loci: CPH4, CPH5, CPH8, CPH12, CPH6, CPH9 [62], CPH22 [63], FH2004, FH2088, FH2096, FH2137, FH2140 [64], CXX.213, C09.250, C20.253 [65] (Table B in S1 File), as their polymorphism was shown in the golden jackal [52]. This number of polymorphic markers is efficient to detect genetic structure and describe genetic diversity within populations [66]. For 12 loci PCR were performed in 15 μl containing 1 μl of DNA, 1 μl of 8 μM primer mix, 7.5 μl of Multiplex PCR (Qiagen). Twelve loci were amplified in three multiplexing sets at following thermal profile: 95°C for 15 min, 40 cycles at 94°C for 30s, 57°C for 90 s, 72°C for 90s and final extension at 72°C for 10 min. The last three loci were amplified individually in total volume of 15 μl containing 1 μl of DNA, 0.5 μl of each 10 μM primers, 7.5 μl PCR Master Mix (EURx). The thermal profile was 95°C for 3 min, 35 cycles at 95°C for 30 s, 57°C for 45 s, 72°C for 45 s and final extension at 72°C for 5 min. PCR products were analyzed in a CEQ8000 sequencer (Beckman Coulter) and allele sizes were estimated using the Beckman Coulter Fragment Analysis Software.

Amplification of hypervariable domain of the mitochondrial DNA (mtDNA) control-region was performed with primers WDLOOPL (5’-TCCCTGACACCCCTACATTC-3’) and H576 (5-CGTTGCGGTCATAGGTGAG-3’) [52]. The PCR reaction mixture containing 2 μl of DNA, 1 μl of each 10 μM primers, 20 μl of PCR Master Mix (EURx) and 16 μl of purified water. The PCR profile was 94°C for 2 min, 40 cycles at 94°C for 15 s, 55°C for 20 s, 72°C for 60 s and final extension at 72°C for 2 min. Amplified products were purified using Clean-up kit (A&A Biotechnology), and then sequenced using BigDye Terminator v3.1 Cycle Sequencing Kit and 3500xL Genetic Analyzer (Applied Biosystems).

Statistical analysis—microsatellites

Polymorphism among microsatellite loci was estimated on three levels. Firstly, we estimated the number of alleles (A), observed heterozygosity (H O), unbiased expected heterozygosity (H E, [67]) and inbreeding coefficient (F IS) for each locus in the total sample (N = 96). The significance of F IS was tested under a randomization procedure, with the Bonferroni correction for multiple comparison. These analyses were performed using GenAlEx version 6.5 [68] and FSTAT version 2.9.3.2 [69]. In addition, a probability test for deviation from the Hardy-Weinberg equilibrium (HWE) was evaluated for each locus using Genepop (Web version 4.2; [70, 71]). Secondly, we estimated polymorphism for each locus in groups of samples designated a priori and corresponding with geographical regions. Aside from A, H O, H E and F IS, we also calculated allelic richness (R; [72]) using FSTAT, as well as mean values for these parameters. HWE was tested for each locus within each region, as well as for each region across all loci. Between-populations genetic differentiation was estimated using F ST [73] as implemented in FSTAT.

To find out whether the geographical grouping of samples corresponded with genetic groups, we applied a Bayesian-clustering method (STRUCTURE version 2.3.4; [74]). Structure was run 15 times for each user-defined number of genetic groups (K = 1–6), with an initial burn-in of 50,000, and 1,000,000 iterations of the total data set. The admixture model of ancestry and the correlated model of allele frequencies were applied. Sampling location was not used as prior information. Next, we examined ΔK statistics that identify the largest change in the estimates of K produced by STRUCTURE (Fig 2A versus Fig 2B) [75]. To visualise the STRUCTURE results we used STRUCTURE HARVESTER 0.6.94 [76]. We then applied CLUMPP 1.1.2 [77] to average the multiple runs given by STRUCTURE and correct for the label switching. The output from CLUMPP was visualised using DISTRUCT v 1.1 [78].

Fig 2

STRUCTURE results: A—estimated likelihoods, ln P(D), of each number of inferred genetic clusters (bars are SD—only given when exceeding the width of dots); B—the corresponding ΔK curves as a function of K; C—ancestry of individuals, estimated for K = 2 and 3 (based on ΔK), and 6 (based on estimated likelihoods).

SEE—south-eastern Europe; CAU—Caucasus; BAL—Baltics; GRE-P—Greece, Peloponnese; GRE-S—Greece, Samos Island.

The Bayesian-based method implemented in the Bayesian Analysis of Population Structure software (BAPS, version 6.0; [79-81]) was used to check the spatial clustering of individuals, and was followed by admixture analysis. In this analysis, geographical coordinates for each sample were used. Ten replicates were run for every upper level of K (2, 3, 4, 5, 10, 15, and 20). The number of iterations used to estimate the admixture coefficient for individuals, and the number of reference individuals from each location were 50 and 200, respectively. The number of iterations applied to estimate the admixture for reference individuals was set at 15.

We also obtained an additional representation of the genetic structure using Principal Components Analysis (PCA). This multivariate descriptive method is not dependent on any model assumption and can thus provide a useful validation of the Bayesian clustering output [82-84]. We used the R package ADEGENET v1.3.4 [85] to carry out the standard PCA. The results of the analysis were presented graphically along first and second axes in line with eigenvalues.

Statistical analysis—mitochondrial DNA

Sequences were aligned in BioEdit software v.7.0.5.3 [86], with alignments then being checked manually. We amplified a 406 base-pair (bp) fragment of the control region for 93 samples also genotyped with microsatellite markers. We were unable to obtain reliable sequences from one sample from Estonia (EST), one from the Caucasus (CAU), and two from south-eastern Europe (SEE). Numbers of haplotypes (H) in the total sample, as well as in particular geographical regions and genetic groups, haplotype diversity (h), nucleotide diversity (π) and mean number of nucleotide differences among haplotypes (k) were all calculated using DNAsp 5.10 [87]. Haplotype frequencies in the overall sample and in each geographical region were calculated using ARLEQUIN v3.5.1.2 [88]. ARLEQUIN was also used to calculate pairwise θ ST among regions using haplotype frequencies. The test for significance was performed with 1,000 permutations. The overall genetic structure, based on haplotype frequencies, was estimated in DNAsp, using H ST ([89]; equation 2). Significance for the global estimate was determined by permutation test, on the basis of 1,000 replicates.

A median-joining haplotype network [90] was constructed in NETWORK v4.6.1.1. (Fluxus Technology Ltd.). We also compared haplotypes identified in this study (GenBank accession nos. KT362174–KT362176) with haplotypes for the golden jackal deposited in GenBank, and originating from Bulgaria, Serbia, Croatia, and Italy (KF588364) [51, 52], Serbia (HQ845260) [91], Bulgaria (AF184048) [92], Poland and Ukraine (KT268318 and KT268319) [93], the Caucasus (KJ490945 and KJ490946) [94], and India (AY289997 and AY289996) [95].

Ethics Statement

Tissue samples used in this study were obtained from individuals that died in vehicle collisions, due to natural causes or as a result of legal hunting. No animal was killed for the purpose of this study.

Results

Microsatellites

From 15 polymorphic microsatellite loci, amplified in 97 golden jackals (Fig 1), we identified 102 alleles (1.05 alleles per individual). At most loci the polymorphism was moderate (5 to 11 alleles). The greatest number of alleles (A = 14) was discovered at locus FH2137, the lowest (A = 3) at CPH5 (Table C in S1 File). In most cases the observed heterozygosity was below 0.60, and at only three loci (FH2004, FH2096, FH2137) did the value exceed 0.70. When all samples were analysed together, 11 of the 15 microsatellites were found not to be in HWE (Table C in S1 File). Similarly, F IS values were found to differ significantly from zero at most loci following Bonferroni correction, the effect being indicative of heterozygote deficiency. Given that all the samples were examined together, and we subsequently found significant substructure, this could be due to the Wahlund effect.

When samples were grouped by geographical distribution, a significant overall F IS was found only in the case of jackals from the Caucasus (CAU). In south-eastern Europe (SEE)—the region represented by the highest number of the samples studied—there were three loci manifesting deviation from the HWE on account of heterozygote deficiency and one, FH2096, indicative of heterozygote excess (Table D in S1 File). SEE was also the only group with significant overall heterozygote deficiency, though F IS was low and non-significant. This group also had the highest mean number of alleles (mean A = 5.40). Allelic richness was similar in SEE, CAU and BAL, though slightly lower in two groups from southern Greece, i.e. from the Peloponnese (GRE-P) and Samos (GRE-S). Observed heterozygosity (H O) was highest in SEE. The lowest H O was found in the insular GRE-S population.

Analysis of genetic structure using Bayesian methods and PCA indicated some grouping patterns. In the STRUCTURE analysis the highest mean likelihood was indicated for six clusters (Fig 2A). GRE-P and GRE-S formed two uniform genetic groups, whereas SEE consisted mainly of individuals from two clusters (with most jackals from Hungary and Romania marked in red, and the majority of those from Serbia and Ukraine shown in violet—Fig 2C; K = 6), but also of individuals of mixed ancestry. Jackals from CAU and BAL were assigned to two other clusters, with more or less equal probability of ancestry from each of them. The ΔK statistic (Fig 2B) suggested two or three genetic groups. In the two-group scenario the first cluster comprised the majority of individuals from SEE and GRE-P, and the second comprised the majority of individuals from GRE-S, CAU and BAL (Fig 2C; K = 2). On the basis of the K = 3 value, BAL and CAU formed the first genetic group, SEE and GRE-P the second, and GRE-S the third (Fig 2C; K = 3). In both of these cases, certain individuals from SEE had the highest probability of ancestry from the CAU/BAL group. These were four individuals from Ukraine (nos. 8599, 8607, 8608, 8927) and two individuals from Serbia (nos. 8620, 8625—Fig 1).

Geographical information about samples in Bayesian analysis (BAPS) suggested the presence of four genetic groups, with a very limited admixture among them (Fig 3A and 3B). In general, the geographical groups designated a priori corresponded to genetic groups as indicated by BAPS. However, one sample from SEE (Ukraine, no. 8608) was assigned to the CAU/BAL cluster, one sample from SEE (northern Greece, no. 8986) was assigned to the GRE-S cluster, and one sample from BAL (Lithuania, no. 9225) was assigned to SEE (Fig 1). A similar result was obtained by way of admixture analysis (Fig 3B), although in this case two additional individuals from SEE (nos. 8927 and 8625 from Ukraine and Serbia, respectively) were found to be of mixed ancestry. Like STRUCTURE, PCA pointed to the genetic distinctness of GRE-S (Fig 4). The remaining samples were divided by PCA into two groups corresponding with SEE/GRE-P and CAU/BAL.

Fig 3

Results of spatial analysis of genetic structure, using BAPS: A—assignment of specimens to four genetic clusters indicated by spatial clustering; B—admixture analysis of identified clusters.

SEE—south-eastern Europe; CAU—Caucasus; BAL—Baltics; GRE-P—Greece, Peloponnese; GRE-S—Greece, Samos Island.

Fig 4

Results of Principal Components Analysis performed in ADEGENET.

First and second axes and corresponding eigenvalues (inset) are shown. SEE—south-eastern Europe; CAU—Caucasus; BAL—Baltics; GRE-P—Greece, Peloponnese; GRE-S—Greece, Samos Island.

Genetic differentiation among the geographical regions was high (overall F ST = 0.199, 95% CI = 0.147–0.258). Pairwise F ST ranged from 0.05 to 0.39 (Table 1). Low genetic differentiation was found between BAL and CAU, whereas all pairwise comparisons with GRE-S indicated a very high level of genetic differentiation (F ST > 0.20). Similarly, marked genetic differentiation was found between GRE-P and CAU and GRE-P and BAL, while moderate genetic differentiation characterised the pairwise comparisons of data for SEE, as set against GRE-P, BAL or CAU.

Table 1

Genetic differentiation among geographical regions: SEE—south-eastern Europe (Croatia, Serbia, Slovenia, Hungary, Romania, Ukraine, northern Greece); CAU—Caucasus (Georgia, Armenia, Mountainous Karabakh); BAL—Baltics (Estonia, Lithuania); GRE-P—Greece, Peloponnese; GRE-S—Greece, Samos Island.

Above diagonal—genetic differentiation calculated from mtDNA haplotype frequencies, below diagonal—genetic differentiation calculated from microsatellites. Significant values (1,000 permutations; P < 0.05) are shown in bold.

Region	SEE	CAU	BAL	GRE-P	GRE-S
SEE		0.347	-0.199	0.507	0.961
CAU	0.125		-0.045	0.024	0.716
BAL	0.100	0.051		0.090	0.863
GRE-P	0.113	0.207	0.268		0.790
GRE-S	0.293	0.205	0.343	0.388

Mitochondrial DNA

Based on the mitochondrial DNA (mtDNA) control region fragment, we identified four unique haplotypes in 93 samples. Both haplotype diversity and nucleotide diversity were low (Table 2), as was the average number of pairwise nucleotide differences (k = 0.706). Apart from BAL, we identified two haplotypes per region. The highest level of haplotype diversity was found in GRE-P, while the most marked nucleotide diversity and highest average number of pairwise nucleotide differences was found in CAU (Table 2).

Table 2

Sample size (N) and genetic characteristics of mtDNA polymorphism in C. aureus in geographical regions and all samples: h—number of identified haplotypes; H [SD]—haplotype diversity and corresponding standard deviation; π [SD]—nucleotide diversity and corresponding standard deviation; k—average number of pairwise nucleotide differences.

SEE—south-eastern Europe; CAU—Caucasus; BAL—Baltics; GRE-P—Greece, Peloponnese; GRE-S—Greece, Samos Island.

Region	N	h	H [SD]	π [SD]	k
SEE	55	2	0.036 [0.035]	0.00009 [0.00009]	0.036
CAU	13	2	0.385 [0.132]	0.00189 [0.00065]	0.769
BAL	4	1	-	-	-
GRE-P	11	2	0.509 [0.101]	0.00125 [0.00025]	0.509
GRE-S	10	2	0.467 [0.132]	0.00115 [0.00032]	0.467
Total	93	4	0.344 [0.061]	0.0017 [0.00033]	0.706

Haplotype H1 proved to be most frequent, being absent only from GRE-S. In BAL this was the only haplotype found. Haplotype H2 proved to be unique to GRE-S, while H3 was shared between CAU and GRE-S, and H4 between SEE and GRE-P (Table 3, Fig 5).

Table 3

Distribution of golden jackal mtDNA haplotypes in the investigated geographical regions.

Frequency in the region and overall frequencies are reported. SEE—south-eastern Europe; CAU—Caucasus; BAL—Baltics; GRE-P—Greece, Peloponnese; GRE-S—Greece, Samos Island.

Haplotype	Motif	SEE	GRE-P	CAU	BAL	GRE-S	Total
H1	TGG	0.98	0.64	0.77	1.00	-	0.800
H2	CAA	-	-	-	-	0.70	0.076
H3	TAA	-	-	0.23	-	0.30	0.068
H4	TAG	0.02	0.36	-	-	-	0.056

Fig 5

The minimum spanning network of mtDNA haplotypes of golden jackals sampled in this study (SEE, CAU, BAL, GRE-P, GRE-S) as well as those desposited in GenBank (Italy, Croatia, Serbia, Bulgaria, Ukraine, NW Poland, the Caucasus, and India).

The length of each line between two circles is proportional to the number of mutations.

Comparing haplotypes identified in this study with those deposited in GenBank (homological sequences of the 250 bp of mtDNA CR), we found that H1 corresponds with a haplotype identified previously in Italy, Croatia, Serbia, Bulgaria, Ukraine, NW Poland, and the Caucasus, while differing by just a single mutation from another haplotype from the Caucasus (H5). A haplotype observed previously in Indian jackals (H6) differed from H4 by just four mutations (Fig 5).

Genetic structure as estimated on the basis of haplotype frequencies was found to be pronounced and significant (H ST = 0.486 for geographical groups, P < 0.001). Pairwise θ ST was highest for the comparison of SEE with GRE-S. No genetic differentiation was noted between BAL and CAU or between BAL and SEE (Table 1).

Discussion

Genetic diversity

Analysing the results obtained with both microsatellite and mitochondrial markers, we found higher genetic diversity than has been reported previously for other European populations of the golden jackal [51, 52], except in the case of the island population from Samos. In Serbia [51] a total of 31 microsatellite alleles at eight loci were found in 120 individuals, giving 3.8 alleles per locus and 0.26 alleles per individual, compared with 6.8 alleles per locus and 1.05 alleles per individual stated in our study (Table C in S1 File). Similarly, in the Serbian population the total observed heterozygosity was 0.28, compared with the 0.52 found in our study. These differences can be explained by the fact that the populations of golden jackals analysed in this study were historically older and larger than those from Serbia, or involved samples from across a larger area, with the SEE geographical group encompassing individuals from a large part of south-eastern, Central and Eastern Europe. Furthermore, the mean number of alleles was higher in SEE and CAU (A = 5.40 and 4.67, respectively; Table C in S1 File) than that found previously [52] in the contemporary samples from Bulgaria (A = 3.5), Slavonia (continental eastern part of Croatia) and Serbia (A = 4.0), Dalmatia (A = 2.8), and Italy (A = 3.7). Moreover, the analysis of the mitochondrial control region revealed four mitochondrial haplotypes (Table 3), as opposed to the one haplotype noted in previous studies [51, 52]. However, the population from Greece (both the Peloponnese and the island of Samos), had a mean number of alleles of around 3.0 (Table C in S1 File), comparable with what was found in Dalmatia [52], and hence slightly lower than the value characterising jackals in Slavonia and Serbia, Bulgaria and Italy [52]. We observed the lowest level of genetic diversity in the island population at Samos (mean A = 2.67, H O = 0.38; Table C in S1 File), which could be explained by the isolation, as low genetic diversity often reflects colonisation of an island by a small number of individuals (the founder effect) and random processes reducing variability, such as genetic drift [96-98].

In the present continent-wide study we supported previous findings of Zachos et al. [51] and Fabbri et al. [52] in Serbia, Bulgaria, Croatia, and Italy, indicating that Europe’s golden jackals harbor less genetic diversity compared to other wild canids, such as wolves [36, 39–41, 99, 100], or red foxes (Vulpes vulpes) [101, 102]. The genetic diversity of European jackals is also clearly lower than that found in jackals from Israel [103], which show signals of hybridization with grey wolves, dogs, and the African golden wolf (Canis anthus) [104]. For example, the mean number of alleles in five populations from Israel ranged from 4.7 to 5.6 versus 2.6–5.4 noted in our study, whereas observed heterozygosity ranged from 0.64 to 0.72 versus 0.38–0.55 in our study. This is despite the dramatic population decline and bottleneck experienced in Israel in the 1960s [55, 105]. Thus the low genetic diversity of Europe’s jackals does not reflect species-specific characteristics, but may be related to the unique history of golden jackals on this continent.

In contrast with the authors of previous studies [51, 52], we noted polymorphisms in the mtDNA control region, even though overall haplotype diversity was low (H = 0.34) with just four haplotypes despite the large sampling area. In jackals from the mainland sites (SEE, CAU, BAL, GRE-P), it was the haplotype recorded in previous studies (H1) that was found to occur most frequently. However, it was absent from the island population (Samos), where the unique H2 haplotype is prevalent. Higher mtDNA diversity compared with previous studies is mainly connected with the larger sampling area including the Caucasus (one ‘new’ haplotype H3) and Greece (three ‘new’ haplotypes: H2, H3 and H4). The highest level of haplotype diversity was found in the Peloponnese. The greater number of haplotypes in the Aegean region could suggest that the present population in Greece may, at least partially, descend from the ancient Greek population. However, confirmation of this hypothesis requires further study, preferably including fossil material. Moreover, in the south-eastern European population, alongside the haplotype discovered by Zachos et al. [51] and Fabbri et al. [52], we identified the additional haplotype H4, which appeared in one animal from the Biruchiy Peninsula (southern Ukraine), i.e., an area outside the Balkan Peninsula, but also in the Peloponnese. Hence, it is possible that the majority of the Balkan population of C. aureus is uniform in regard to control-region polymorphisms, as suggested by earlier studies [51, 52]. Hence, despite the discovery of additional haplotypes, the genetic diversity in the mitochondrial control region in Europe’s golden jackals should be regarded as low when compared with that in other canids [102, 106, 107]. However, further sampling will probably result in the detection of new polymorphisms in mtDNA of the golden jackal as a species, as the haplotype found in Indian jackals (denoted as H6 in Fig 5) differs by 4–6 substitutions from the haplotypes identified in the present study.

Genetic structure

Previous studies of golden jackals in Europe emphasized the limited degree of genetic structuring, with only the coastal population from Dalmatia clearly differentiated from other Balkan samples [51, 52]. A genetic identity relating to Dalmatia has also been suggested in the case of the grey wolf [108], and was explained either by reference to an origin of this population in a distinct refugium, or in terms of ecological and behavioural factors [41, 109, 110]. Fabbri et al. [52] also noted markedly smaller number of alleles (A = 2.8) and more limited heterozygosity (H O = 0.37) in Dalmatian jackals and suggested a long-term isolation of this population. In respect to this, we also call to attention that golden jackals were present in southern Dalmatia already in the Middle Ages [111] and possibly even much earlier [58].

Our analysis extending to the whole of Europe has pointed to the existence of a pronounced genetic structure in relation to both nuclear and mitochondrial markers. Individuals from an extensive area of south-eastern Europe generally form a uniform genetic group, as already noted by Zachos et al. [51]. Fabbri et al. [52] also reported small genetic differentiation in microsatellite markers among populations from Bulgaria, Slavonia and Serbia. This probably reflects recent expansion of the species in this region. However, Greek samples indicate the existence of a distinct population in the Peloponnese (STRUCTURE, BAPS) (see also [112]), even if both haplotypes found in south-eastern Europe were also present in animals from this peninsula. We can speculate that our results support the hypothesis that an ancient Greek population survived in the Peloponnese to the present day, recently merging with a population expanding in from the east. A similar interpretation can be put forward in regard to Dalmatian jackals, as already suggested by Fabbri et al. [52]. Thus the two known areas with the early Holocene findings of jackals [58, 59] are also the only two areas in south-eastern Europe today that show higher genetic differentiation, giving further support for the continuous presence of ancient populations along the Mediterranean coast.

STRUCTURE and BAPS suggested ongoing gene flow between the Caucasus and Europe as well—some individuals from SEE had the highest probability of ancestry from the CAU/BAL cluster. Interestingly, when the microsatellite genotypes are concerned, an individual from south-eastern Europe (no. 8927; Table A in S1 File) with the additional haplotype H4 (which is frequently found in the Peloponnesus Peninsula), was identified as having ancestry from the Caucasus (STRUCTURE: two-clade and three-clade scenarios) or mixed ancestry from the Caucasus and Samos Island (BAPS: admixture analysis).

The island population of golden jackals on Samos was highly differentiated from those from other sampling sites (F ST, STRUCTURE: three-clade scenario, PCA, BAPS). Unfortunately, there was no access to samples from the Turkish mainland, so it remains unclear whether the geographical barrier of water restrics gene flow between the island and that mainland. However, our genetic data indicate that there are or were some connections between Samos and northern Greece (e.g. an individual no. 8986 sampled in Chalkidiki Peninsula [Fig 1, Table A in S1 File] was assigned to the Samos cluster).

Stepping stone model or long-distance colonizers—on the origin of the Baltic jackals

First possible observations of jackals in the Baltics are known from 2011, when groups of several jackals were noted in Estonia [113]. In 2013 and 2014 several animals were shot, photographed, or detected during howling surveys in Estonia and Latvia [113, 114], and in 2015 the first jackal was shot in Lithuania [115]. Although several carnivore experts suggested that natural expansion was likely, the governments of the Baltic States decided to assume that jackals were introduced by humans [38].

The genetic data suggest that jackals from the Baltics originate from the Caucasus region (Estonian samples), and from the population expanding out of south-eastern Europe (Lithuanian case). This dual origin does not support the idea that jackals were introduced by humans, as it is unlikely that someone would capture jackals in different regions and smuggle them to the Baltics. Additionally, recent records of jackal occurrence from Slovakia, Ukraine, Belarus, and north-western and eastern Poland [38, 93, 113–116], suggest that both Caucasian and southeastern European populations are spreading towards the north. The presence of the Caucasian gene pool was also detected in animals from NE Ukraine, further supporting the hypothesis of natural expansion from the Caucasian region through Ukraine towards Estonia.

The dynamics of species’ range expansions depend on habitat connectivity, but also on dispersal ability [118] and habitat plasticity [119]. Two basic models were suggested for a dispersal through fragmented environment, where suitable habitat is distributed as a series of patches. In the ‘island model’ all patches are equally accessible, while in the stepping-stone model exchanges of individuals are restricted to adjacent populations [120]. Although previous genetic data suggested a ‘stepping-stone’ nature of golden jackal dispersal [52], our results indicate the possibility of long-distance dispersal in this species. This can also be supported by a review of literature data, which includes several records of sudden appearances of jackals far from other known populations. Such example include the (re)colonization of Hungary in the 19th century [60] with the closest known populations at that time being in Dalmatia, Croatia (at ca. a 300 km straight-line distance) or Bulgaria (400 km away).

Another case resembling the sudden occurrence of several jackals in the Baltics, refers to the first colonization of Slovenia in the mid-20th century. In winter 1952/3 several jackals suddenly appeared in Central Slovenia near Ljubljana, with reported observations of groups of up to six animals [121] and later shooting of two animals near Ljubljana and one at the foothills of the Julian Alps in NW Slovenia [122]. At that time, the closest jackal population was known from Ravni Kotari in Dalmatia, Croatia [123], approximately 210 km from Ljubljana. In this probable case of long-distance dispersal, jackals seem to have dispersed in a group, as it would be highly unlikely that several animals would appear independently at the same time in the same place so far from the closest population.

More recent records that can be considered potential cases of long-distance dispersal of jackals include:

a male observed several times from 1996 and then shot in 1998 in Südbrandenburg in Germany [124] and an individual photographed four times in 2012 in the Bavarian Forest [125]. These records were 430 km and 270 km distant, respectively, from the closest-known reproducing population in eastern Austria and western Hungary;

five photo-records of a jackal in 2011 in the Northwestern Alps of Switzerland [126], with the closest known reproducing population in NE Italy 450 km away;

an individual shot in 2014 near Olevsk in Northern Ukraine [127], 430 km from the closest known population in Southern Ukraine and Moldova;

an individual shot in 2012 near Tomašovka in Belarus [117], 410 km from the closest known reproducing population in Hungary;

a young male found dead in April 2015 on a road in NW Poland, close to the German border [98], ca. 610 km from the reproducing population in NW Hungary;

GPS-GSM collared 1.5 year old female, which travelled 220 km during 12 days in Hungary in 2014 (J. Lanszki unpubl. data).

Based on this review of jackal occurrences and our genetic data, we suggest that it is not uncommon for golden jackals to disperse over several hundred kilometers in human-dominated landscapes. This could explain the speed of jackal expansion in Europe that has been observed in the last decades [38]. We also suggest that the recent colonization of the Baltic States is most likely a case of long-distance dispersal. The first ‘wave’ of colonization of the Baltics appears to have originated from the Caucasus region via Ukraine. The second wave on the other hand seems to have originated from south-eastern Europe through an expansion front in Romania, Hungary/western Ukraine, Slovakia, and Poland. According to available records it even appears that a group of several jackals can disperse together (see also [116]). If true, this would have important implications, as it would considerably increase probability of successful colonization of new areas.

Management and conservation implications

The golden jackal has already been declared an alien, potentially invasive species in all Baltic States (e.g. [128]). However, an Invasive Alien Species (IAS) needs to meet at least three criteria: 1) it should be non-native, allochtonous, introduced by people; 2) it should threaten biological diversity on the local scale; and 3) it should be characterised by rapid population growth [129]. Although exponential increase in population size has been observed (e.g. in Hungary [130]), the other two criteria have not been met. The movement north is evidently a result of natural migration (as the present study shows), and there is no proof of a harmful effect on local fauna [131-134]. Also there are no major complaints about golden jackals inflicting harm on domestic animals reported from Europe [133-136]. Occasional reported claims of jackal depredation of livestock are believed to be exaggerated often [134, 137], or connected with erroneous identification, when reported cases have been inspected using forensic genetics [138]. Recent genetic analysis [104] has also shown that the severe impacts on livestock reported from Israel [139], are probably not connected with golden jackals per se, but rather with individuals of admixed origin between several canid species. Furthemore, the parasite load in the European golden jackal is similar to or lower than that in other carnivores (e.g. the red fox, grey wolf, and wild cat [Felis silvestris]) in the region [140-143], and no attacks by jackals on people are known. For these various reasons, concerns regarding serious negative impacts of the expansion of the golden jackal in Europe appear to be unfounded as yet.

Nevertheless, results presented here have several management and conservation implications. The existence of long-distance dispersal in the golden jackal would seem to warrant the initiation of international coordination in management of the species in Europe and more focus on management at the population, rather than at the national level, especially considering considerable differences that currently exist among countries [38]. We therefore suggest the development of trans-boundary management strategies and documents similar to the population-level management approaches developed in the case of Europe’s large carnivores [144]. We also call for a revision of the approach used in managing jackals in the Baltic States, given that our results contradict the presumption of the local decision-makers about the human-assisted origin of the Baltic population. Lastly, our results provide a basis for the development of a conservation strategy for the golden jackal in the region. We propose that priority should be given to the Caucasus region, which harbors high genetic diversity in terms of the number of microsatellite alleles, as well as to the regions of the Peloponnese and Dalmatia [52], in which a relict gene pool from ancient Mediterranean populations appears to have persisted. The golden jackal is listed as an Annex V species in the EU Habitats Directive and as such, taking above into account, should be legally protected in all EU member states (for legal implications of range expansion in this species see [38]).

Material studied and its genetic characteristics.

Table A, List of examined specimens including specimen number, sex, locality information, date, geographical coordinates, and mtDNA haplotype designation. Table B, Microsatellite genotypes. Table C, Per locus genetic diversity in 97 samples of golden jackal. Table D, Per region genetic diversity estimated based on polymorphisms in 15 microsatellite loci.

(DOC)

Click here for additional data file.