Using genetic tools to estimate the prevalence of non‐native red deer (Cervus elaphus) in a Western European population

Abstract Game species like the red deer have been subjected to anthropogenic impacts for centuries. Translocations are often carried out—sometimes illegally—not only for sporting purposes, but also to increase trophy quality, reduce inbreeding, or mitigate bottlenecks after excessive persecution. Apart from the blurring of large‐scale genetic structure, translocations without adequate quarantine measure risk introducing pathogens into potentially immunologically naïve populations. It is therefore important to understand the frequency of clandestine translocations. Identification of non‐autochthonous animals and their potential origin is often difficult and, in red deer, has been hampered by the lack of large‐scale genotypic datasets for comparison. In the present study, we make use of a recently published European‐wide microsatellite dataset to detect and quantify the presence of non‐autochthonous red deer in a large population sample (n = 1,780) from Central Europe (Belgium). Using factorial correspondence analysis, assignment tests and Bayesian clustering algorithms we arrive at an estimate of 3.7% non‐autochthonous animals (or their descendants). Some of these animals were assigned to a nearby French population and may have immigrated into Belgium naturally, but the large majority must have been introduced by humans. Our analysis pointed to the British Isles and Germany/Poland as the potential origin of many introduced deer, regions known to have been source populations for translocations in Europe and beyond. We found evidence for recreational hunters using carcasses from farmed deer to fulfill mandatory hunting quotas. Our study is the first to quantify the extent of human‐mediated introductions in a European game species at such a large scale with large and representative sample sizes.

it remains unclear for many regions whether the occurrence of both lineages is natural or due to reintroduction of red deer after local extirpation. The large-scale pattern of the three lineages that are geographically separated, however, strongly suggests that the natural genetic structure at the level of glacial refugial lineages is still present in red deer.
To date, there are few studies directly aiming at identifying the genetic signature of translocation events. Probably, the most extensive such study was that of Frantz et al. (2006) who were able to prove the illegal introduction of a small number of deer into a local population in Luxembourg. In general, this kind of analysis is difficult as human-mediated introductions cannot be unequivocally distinguished from natural immigration from a neighboring, but genetically differentiated population. Also, if the origin of a translocated animal is to be determined, DNA samples from ideally all possible source populations need to be available. Because until recently no European-wide dataset based on high-resolution markers such as microsatellites was available, genetic studies on translocations have long been hampered. In this study, we present an in-depth analysis of regional red deer populations in Belgium-where since 1994 no private individual is legally allowed to translocate wildlife-based on 13 microsatellite loci and 1780 samples and compare this dataset to the recently published European-wide dataset (Zachos et al., 2016) to identify nonindigenous genotypes and to arrive at a quantitative estimate on how large the proportion of potentially translocated individuals is, thus testing claims that few if any European red deer populations are free from introductions (Hartl et al., 2003).

| Sample collection and laboratory work
Between 2003 and 2009, we collected tissue samples from 1,733 red deer from the southern Walloon part of Belgium, as well as 47 samples from the neighboring Eifel region of Germany ( Fig. 1; see also Frantz et al., 2012 who used about half of our present samples). For ease of reference, we will refer to the dataset as "Belgian" or "Walloon," despite the few German animals. The total forested area of Wallonia amounts to approximately 5,000 km 2 , which corresponds to one-third of the total area of the region. Human population density in the southern part of Wallonia, where our study was based, was mostly between 20 and 90 inhabitants⁄km 2 , depending on the municipality (Thomas, Frankhauser, & Biernacki, 2008). In 2014, the Walloon red deer population was estimated to have contained approximately 11,000 individuals (Direction de l'Etat Environnemental 2014). Samples were collected during legal hunts from harvested animals. The center of the forest management unit (a so-called "triage") where an individual was harvested was recorded as its sampling location.
DNA was extracted using a chloroform-based extraction method (Doyle & Doyle, 1990). Samples were genotyped using 13 microsatellite loci (BM1818, Cer14, CSPS115, CSSM14, CSSM16, CSSM19,   CSSM22, CSSM66, ETH222, Haut14, ILSTS06, INRA35, and MM12; for references see Kuehn et al., 2003) in three multiplex polymerase chain reactions (PCR) using the Qiagen Multiplex kit (Qiagen, Hilden, Germany). Detailed information on the PCR composition and reaction times can be found in Dellicour et al. (2011). We obtained a complete 13-loci genetic profile for 1,746 of the 1,780 deer, with the remaining profiles having genotypes at a minimum of 11 loci. In all, 13 animals had a missing genotype at locus ILSTS06. Reactions were performed using a Verity thermocycler (Applied Biosystems, Warrington, UK). PCR products were separated using an ABI 3100 automated DNA sequencer (Applied Biosystems), and the data were analyzed using GeneMapper version 3.7 (Applied Biosystems).

| European reference data
We aimed to identify the source population of translocated deer, both to identify illegal introductions and animals that had dispersed naturally, using a European microsatellite reference dataset analyzed by Zachos et al. (2016). Not considering Belgian animals, it contained genetic profiles of 608 red deer from 26 locations throughout the continent (see Fig. 2), including 30 samples from a French deer farm (Boisgervilly). Both the Walloon and the European dataset were genotyped in the same laboratory using the same markers. We obtained a complete 13-loci genetic profile for 571 of the 608 European deer, with the remaining profiles having genotypes at a minimum of 11 loci. Loci Cer14 and ILSTS06 did not amplify in 14 and 10 individuals, respectively. Using the individual-based modal population mixture analysis implemented in BAPS v.5.4 (Corander, Waldmann, Marttinen, & Sillanpää, 2004), the European dataset was inferred to consist of 25 genetic populations (Zachos et al., 2016). Seven of these 25 clusters consisted of six individuals or less. These were excluded from the present analysis. Although the animals sampled at a deer farm were assigned to different genetic clusters, for the purposes of the present study, we created a separate reference population containing all the farmed deer. The European reference dataset therefore consisted of 596 non-Belgian individuals belonging to 19 genetic clusters (18 BAPS-defined partitions with N > 6, plus the deer farm; Fig. 2).

| Translocation analysis of Belgian red deer
While not having any a priori suspects, we needed to remove genetic outliers to avoid their interference in the analysis of the population genetic structure of Wallonia, as well as in the calculation of exclusion probabilities. We therefore used GENETIX v.4.05.2 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 2004) to perform a factorial correspondence analysis (FCA) to visualize the genetic distance between all 1780 tested red deer. The coordinates of the two-first principal components were entered into ARCVIEW 3.3 (ESRI Inc., Redlands, CA, F I G U R E 1 Geographic origin of Belgian red deer samples included in this study and location of the genetic subpopulations inferred using the STRUCTURE (top) and spatial BAPS (bottom) algorithms. The size of the pie charts indicates the number of samples collected from a locality, whereas the pattern of the pie chart indicates the identity of the genetic clusters USA) and we used the Animal Movement Analyst Extension (Hooge & Eichenlaub, 1997) to remove the outlying 5% (N = 89) of all individuals using the harmonic mean method (Dixon & Chapman, 1980). The individuals that were outliers in the FCA analysis were excluded from the analysis of the population genetic structure of red deer in Belgium (it was later formally tested whether these deer were nonautochthonous using an exclusion test; see below).
We analyzed the population genetic structure of the whole Belgian dataset (minus outliers) using two different (Bayesian) genetic clustering algorithms. First, we analyzed the data using the program STRUCTURE v. 2.3.1 (Pritchard, Stephens, & Donnelly, 2000).
To estimate the number of subpopulations (K), 10 independent runs of K = 1-10 were carried out with 10 6 Markov chain Monte Carlo (MCMC) iterations after a burn-in period of 10 5 iterations, using the model with correlated allele frequencies and assuming admixture.
ALPHA, the Dirichlet parameter for the degree of admixture, was allowed to vary between clusters. After deciding on the most probable number of sub-populations based on the log-likelihood values (and their convergence) associated with each K, we calculated each individual's percentage of membership (q) for each cluster, averaging q over 10 runs. Individuals were assigned to the STRUCTURE cluster for which they had the highest q value. Finally, we also analyzed the data using the spatially explicit genetic clustering method that is implemented in the program BAPS v.6.0 (Corander, Sirén, & Arjas, 2008).
In addition to the genetic data, the algorithm considers the specific geographic coordinates of each individual and modally assigns each individual to its population of origin. We performed 10 runs for each K = 2-10. We calculated the average assignment proportions for each sampling location and mapped them using ArcGIS 10.3 (ESRI Inc., Redlands, CA, USA).
For each STRUCTURE-and BAPS-defined cluster, we tested for the significance of heterozygote deficiency or excess using the Markov chain method in GENEPOP v.4.0 (Raymond & Rousset, 1995), with 10,000 dememorization steps, 500 batches and 10,000 subsequent iterations. The false discovery rate technique was used to eliminate false assignment of significance by chance (Verhoeven, Simonsen, & McIntyre, 2005). We tested each STRUCTURE cluster for linkage disequilibria among loci using an exact test based on a Markov chain method as implemented in GENEPOP 3.4. Since the presence of immigrants can lead to linkage disequilibria among loci can (Paetkau, Slade, Burden, & Estoup, 2004), we did not exclude linked loci from further analysis.
Throughout this study, we used GENECLASS 2.0.g (Piry et al., 2004) to calculate the probability of an animal belonging to a genetic population (exclusion probability) based on the Monte Carlo method of Paetkau et al. (2004). We simulated 10,000 multi-locus genotypes and set the threshold for individual exclusion to 0.01. In wildlife forensics, a more stringent threshold for excluding animals from a population-such as p < .001-is considered necessary (Manel, Berthier, & Luikart, 2002), but an exclusion threshold of p < .01 is normally used in ecological studies to identify genetic immigrants (e.g., Aspi, Roininen, Ruokonen, Kojola, & Vila, 2006;Clark, Brown, Stechert, & Zamudio, 2008;Proctor, McLellan, & Barclay, 2005). Individuals were assigned to their most likely source population (assignment test) using the partial F I G U R E 2 Geographic location of the European red deer reference populations and composition of the genetic populations inferred using the individual-based BAPS algorithm from Zachos et al. (2016). The size of the pie charts indicates the number of samples collected from a locality, while the pattern of the pie chart indicates the identity of the genetic clusters. "Mesola II" was excluded from the analysis as it contained only six individuals. The entire deer farm (not indicative of geographic location) was considered to be a distinct reference population although its animals were assigned to different BAPS clusters. The locations of the Walloon red deer are indicated by individual sampling locations Bayesian approach of Rannala and Mountain (1997) implemented in GENECLASS.
In order to avoid false inclusion of non-autochthonous individuals in the Belgian reference dataset (not including the 89 outliers), we first calculated the exclusion probability of each Walloon deer for the three STRUCTURE-defined reference populations using a leaveone-out approach where each individual in turn is excluded from the populations during computation. Furthermore, for each Belgian deer (that was not an outlier), we calculated the probability of belonging to each European reference population and assigned it to its most likely Belgian or European source population. The following Belgian animals were considered to be non-autochthonous and excluded from the reference dataset for further analysis: (1) deer that could be excluded with p < .01 from all STRUCTURE-defined Belgian reference clusters, (2) animals that were assigned with a (total) score of >0.85 to one (or two) European reference population(s). Furthermore, reference animals that were assigned with a score of <0.85 to a Belgian population and with a score of >0.15 to a European reference population were considered to have at least one non-autochthonous recent ancestor.
We choose the 0.85/0.15-threshold as a compromise between correctly estimating the number of admixed individuals (efficiency) and identifying admixed individuals (accuracy). According to simulations by Vähä and Primmer (2006), the accuracy of detecting admixed individuals increases with increasing q value, but the total number of admixed individuals is over-estimated. Given the number of loci used (N = 13) and the degree of genetic differentiation between the reference populations (see Zachos et al., 2016), a 0.90/0.10-threshold would be most likely to provide a correct estimate for the number of admixed individuals, whereas a 0.80/0.20-threshold would allow identification of all admixed individuals (Vähä & Primmer, 2006).
Finally, considering the same criteria to assess non-native status, we identified non-autochthonous individuals among the 89 FCAoutliers by first calculating the exclusion probability of each animal for every reference cluster (Belgium and Europe) and then assigning them to their most likely population of origin. The harvest locations of the non-native deer were mapped using ArcGIS 10.3.

| RESULTS
After excluding the 89 outliers identified using an FCA (Fig. 3), the loglikelihood values generated by the STRUCTURE analysis supported the presence of three genetic clusters in the Belgian dataset (Fig. S1,   Fig. 1). Similarly, the individual-based spatial model implemented in BAPS inferred the presence of three geographically coherent genetic clusters in red deer from Wallonia (Fig. 1). Both programs roughly infer the same clustering solution. The spatial BAPS algorithm identified the deer sampled in the eastern half of the study area to form one genetic cluster. The deer in the western half formed two distinct genetic populations with the genetic discontinuities corresponding to a major motorway bisecting the study area (Fig. 1). The non-spatial STRUCTURE algorithm suggested a more gradual differentiation between the three clusters.
While one to two loci deviated from Hardy-Weinberg equilibrium (HWE) in each STRUCTURE-defined cluster after correcting for multiple tests, this was the case with zero to three loci in case of the BAPS-derived clusters (Table S1). Locus CSSM66 deviated strongly from HWE in two of the three clusters, independently of the algorithm considered. It was therefore excluded from further analysis, which were based on the STRUCTURE clusters, since fewer loci overall deviated from HWE when considering the partition results of this algorithm (Table S1). After correcting for multiple tests, between one and three pairs of loci were in linkage disequilibrium in each STRUCTUREdefined cluster.
Using the leave-one-out exclusion approach, eight deer from the Belgian reference dataset (89 outliers not included) could be excluded (p < .01) from all three STRUCTURE-defined reference populations.
Of these, seven could be excluded from all 19 European reference populations (Table 1; suggesting they originated from a non-sampled population), while the eighth animal could be assigned to the NE German/E Poland reference cluster. A further five individuals could not be excluded from the Belgian reference clusters but were assigned to one (or two) European reference populations. Two of these five were assigned to the Eastern France population and could possibly have migrated naturally into the study area. Finally, 10 Belgian reference individuals were shown to be admixed with another European reference population (Table 1). Three of these 10 were admixed with Eastern France, the nearest European reference population. All these 23 (8 + 5+10) genotypes were excluded from the Belgian reference populations for further analysis.
Altogether 40 of the 89 FCA-outliers could be excluded from all three Walloon clusters at the p < .01-level (28 animals even with F I G U R E 3 Factorial correspondence analysis of Belgian red deer (N = 1,780). The analysis was based on 13 microsatellite loci. The 89 (5%) outliers were identified using a harmonic mean method. The percentage of the total variation explained by each of the two axes is given. The inset magnifies the interface between the core population and the outliers T A B L E 1 Identification of non-autochthonous Belgian deer. We first removed 89 deer from the Belgian dataset that were outliers (5%) in an FCA analysis. For each animal, we then calculated the probabilities of it belonging to each of the three Belgian reference clusters dataset and to each of the 19 European reference populations by means of assignment tests with the GENECLASS software. Animals that could not be excluded from all reference populations at the p < .01-level were assigned to their most likely population of origin. We considered the two mostly likely populations of origin if the assignment score for the first population was <85%. For further details, see Section 2. (a) Non-outlying Belgian deer that could be excluded from the three Belgian populations using a leave-one-out approach; (b) non-outlying Belgian deer that were assigned with confidence to a European reference population; (c) non-outlying Belgian deer that had a recent non-native ancestor; (d) Belgian outliers that could be excluded from the three Belgian reference populations; (e) non-excluded Belgian outliers that could be assigned with confidence to a European reference population. Max. Belgium = the maximum exclusion probability observed in any of the three Belgian reference clusters. Max. Exclusion Europe = the maximum exclusion probability observed in any of the 19 European reference clusters. Animals in bold may have migrated naturally into the study area. The geographic location of the reference clusters is given in Fig. 2  therefore identified as being non-autochthonous or as having recent non-autochthonous ancestors. A small population of deer in the northeast of the study area some distance away from the nearest regular population appears to be entirely non-autochthonous (Fig. 4). One non-native individual was observed in the German part of our study region. Generally, the non-native animals occurred across most of the

ID
study area (few non-autochthonous deer were observed in southwestern Wallonia) rather than being concentrated in a few localities (Fig. 4).
While non-autochthonous deer were harvested during the whole hunting season (21st Sep. to 31st Dec.), there was a clear increase in non-natives toward closed season (Fig. 5). This increase was not observed with native deer.

| DISCUSSION
Translocations of game animals are often carried out, sometimes illegally (see Frantz et al., 2006), for sporting purposes and to increase trophy quality (in the case of red deer: body mass, size, and number of points in antlers), but apart from blurring natural genetic patterns the potential transmission of diseases is also an imminent danger (Ferroglio et al., 2011). Until recently, the unequivocal identification of non-autochthonous red deer and their potential origin within a geographically limited region was impeded by the lack of a continent-wide nuclear analysis of the species and difficulties in inter-laboratory comparisons of microsatellite alleles. We genotyped 1,780 red deer from Belgium and the Eifel at 13 microsatellite loci and compared these profiles with a European red deer dataset. The combined data of almost 2,500 red deer microsatellite profiles enabled us for the first time to quantitatively assess the proportion of non-autochthonous red deer within a geographically limited region based on nuclear genetic markers. Given the absence of a European reference dataset, Frantz et al. (2006) were not able to confidently identify anthropogenic introductions. The present study is therefore, to our knowledge, the first to address this topic in such detail and with such an extensive dataset.
Our dataset confirmed the conclusions by Frantz et al. (2012) that the E411 motorway in the west of our study area represented a gene flow barrier to red deer. The motorways in the eastern part of our Belgium study area did not appear to have the same effect (Fig. 2).
A more detailed landscape genetic study will be necessary to further untangle the effect of geographic distance and landscape features on the population genetic structure of Belgian red deer. The information provided in this study will help to avoid analytical problems related to the inclusion of non-autochthonous animals.
The inferred frequency of non-autochthonous individuals was between three and four percent. While this value seems rather low, it corresponds roughly to one in 25-30 red deer, which in total amounts to a large number of animals-approx. 400 based on a census size of roughly 11,000 head in Wallonia. Some of these non-native animals may have immigrated naturally from adjoining populations (e.g., F I G U R E 4 Geographic location of nonnative deer (N = 66). Red: animals excluded from all three Walloon STRUCTURE clusters (see Fig. 1 We may have underestimated translocation activities, since our marker system is unlikely to have enough resolution to identify all non-natives or descendants of non-native animals. The markers per se have proved suitable for the study of both small-scale and large-scale genetic structuring and for the identification of immigration/introduction events (Dellicour et al., 2011;Frantz et al., 2008Frantz et al., , 2012Zachos et al., 2016). Still, analyzing a higher number of loci may have pushed a (small) number of individuals below an exclusion threshold. More generally, a microsatellite dataset will only yield information on events involving a few recent generations, while incidents further back in time are unlikely to be detected (e.g., Frantz et al., 2013). It is known that red deer have been translocated for many centuries (Apollonio et al., 2014;Niethammer, 1963), and in many cases these interferences with the natural distribution pattern and its demographic, genetic, and evolutionary consequences cannot be uncovered anymore by means of molecular (or any other) approaches.
It is not clear to what extent our thresholds (FCA-outliers, exclusion, and assignment tests) for identifying non-autochthonous animals were appropriate to produce an accurate estimate of the number of non-native deer. We failed to detect some non-autochthonous animals by excluding the 5%-FCA-outliers from the dataset. However, our results suggest that, in principle, the approach allows the identification of introduced animals in the absence of a priori suspects (which was the case in Frantz et al., 2006). Removing a larger percentage of all individuals might be necessary to ensure that no non-native animals remain non-included. Only considering animals with an exclusion probability of p < .001 as being non-native would have reduced our estimate of the number of non-natives, but would in all likelihood have been too conservative (e.g., Frantz et al., 2006). Considering animals that were assigned with a score of >0.85 to a Belgian population appeared a good compromise between accuracy and efficiency in the identification of admixed individuals, but may have overestimated the number of admixed individuals while not identifying all admixed individuals (Vähä & Primmer, 2006). However, while an individual threshold might be subject to debate, most animals fulfilled more than one criterion that classified them as non-native.
Our analyses revealed two populations as the most likely source Scottish deer on the one hand and deer from NE Germany/E Poland on the other (see Fig. 2). As a result, based on the available data, direct introductions from deer farms-and therefore indirectly from the  (Cunningham, 1996;Woodford & Rossiter, 1993). While in our case it might not be possible to differentiate between intentional introduction and unintentional escape of red deer from Walloon deer farms, both are likely to increase the risk of disease introduction.
If our results were representative for other regions-at least with respect to the order of magnitude-they would indicate a substantial anthropogenic impact on populations of one of the most widespread large European mammals. While results by Frantz et al. (2006) suggest that the results in the present study might be high compared to other regions, a few south-eastern populations in our European reference dataset appeared to contain non-autochthonous animals (see  (Ludt et al., 2004;Niedziałkowska et al., 2011;Skog et al., 2009). The bi-parentally inherited nuclear microsatellites were expected to show less pronounced geographic structuring across Europe because of their higher mutation rates that result in a faster erosion of postglacial demographic signals than with mtDNA.
Surprisingly, Zachos et al. (2016) found a clear signal of three genetic clusters at the highest hierarchical level in their multi-locus dataset that was geographically congruent with the one known from mitochondrial DNA phylogeography.
Our present study with high enough resolution power (i.e., a sufficient number of loci and large sample sizes) suggests that, at smaller scales than those usually addressed in phylogeographic analyses, natural structures have been blurred to a certain extent by the introduction of non-native individuals that subsequently reproduced successfully.
This is likely to hold true not just for the Belgian red deer, but also for other European populations, including those of other ungulate species.

ACKNOWLEDGMENTS
This work was supported by grants from the Public Service of Wallonia

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTIONS
ACF, MCF, and SB designed the study. SB and MC collected samples.
MCE and MCF generated the raw data in the laboratory. ACF, FEZ, and MCF analyzed the data. All authors contributed to the writing of the manuscript.