Genetic Structure and Gene Flow in Red Foxes (Vulpes vulpes) in Scandinavia: Implications for the Potential Future Spread of Echinococcus multilocularis Tapeworm

Knowledge about the dispersal and gene flow patterns in wild animals are important for our understanding of population ecology and the connectedness of populations. It is also important for management relating to disease control and the transmission of new and emerging diseases. Our study aimed to evaluate the genetic structuring among comparative samples of red foxes in a small part of Scandinavia and to estimate the gene flow and potential directionality in the movements of foxes using an optimized set of microsatellite markers. We compared genetic samples of red foxes (Vulpes vulpes) from two areas in Sweden and two areas in Norway, including red fox samples from areas where the occurrence of the cyclophyllic tapeworm Echinococcus multilocularis has been documented, and areas without known occurrence of the parasite. Our results show a high level of gene flow over considerable distances and substantiates migration from areas affected with E. multilocularis into Norway where the parasite is not yet detected. The results allow us to better understand the gene flow and directionality in the movement patterns of red foxes, which is important for wildlife management authorities regarding the spread of E. multilocularis.


Introduction
Dispersal is a key element that connects animals and populations through gene flow [1]. However, the importance of the role of dispersal in connecting populations is still not fully understood [2].
A total of 150 individual red fox samples were used in this study, with 101 sampled from Hedmark County in Norway and 49 sampled from Dalarna, Södermanland, and Östergötland Counties in Sweden (Table 1, Figure 1). Individual red fox samples were grouped into the nearest geographical area to ensure a large enough sample size. Samples from Hedmark South consisted of 52 individuals and samples from Hedmark North consisted of 49 individuals. For Sweden, 24 samples were collected for Sweden North and 25 samples were collected for Sweden South. Red fox samples were either collected as hair from live trapped individuals (n = 33) or collected as tissue samples from shot foxes (n = 117). The sampling was opportunistic in that we collected material from ongoing studies. Sampling was performed from 2012 to 2013.
DNA was extracted from selected hair root or dissected muscle samples according to standard protocols [28] involving lysis with proteinase K, sodium dodecyl sulfate (SDS), phenol/chloroform extraction, and ethanol precipitation. Other dissected muscle samples were processed for DNA A total of 150 individual red fox samples were used in this study, with 101 sampled from Hedmark County in Norway and 49 sampled from Dalarna, Södermanland, and Östergötland Counties in Sweden (Table 1, Figure 1). Individual red fox samples were grouped into the nearest geographical area to ensure a large enough sample size. Samples from Hedmark South consisted of 52 individuals and samples from Hedmark North consisted of 49 individuals. For Sweden, 24 samples were collected for Sweden North and 25 samples were collected for Sweden South. Red fox samples were either collected as hair from live trapped individuals (n = 33) or collected as tissue samples from shot foxes (n = 117). The sampling was opportunistic in that we collected material from ongoing studies. Sampling was performed from 2012 to 2013.
DNA was extracted from selected hair root or dissected muscle samples according to standard protocols [28] involving lysis with proteinase K, sodium dodecyl sulfate (SDS), phenol/chloroform extraction, and ethanol precipitation. Other dissected muscle samples were processed for DNA extraction according to the manufacturer's instructions using the MasterPure Complete DNA and RNA Purification Kit (Epicentre, Madison, WI, USA). Table 1. Sampling information for the 150 red fox individuals from Norway and Sweden included in this study. Columns describes study area, country, total fox analyzed and separated into sex, the sampling method (GPS fox-trapped and released alive; shot fox). See also Figure 1 for geographical localization of all individual samples. Fifteen established microsatellite markers and one polymorphic sex-linked locus previously employed in multiplex [29][30][31][32][33] and 11 relatively new markers tested before only in uniplex [33] were processed groupwise using the MPprimer software package [34] to design primers for increased degrees of multiplexing. Newly designed primers for all 27 loci were tested and optimized for use in three multiplex panels consisting of 11, 10, and 6 primer sets. Marker and primer information is provided in Supplementary Table S1. All multiplex PCR reactions using the final primer concentrations listed in Table S1 were composed as follows: 1 × Buffer B1, 1.5 mM MgCl 2 , 0.1 mM deoxyribonucleoside triphosphates (dNTPs), 0.1 mg/mL bovine serum albumin (BSA), 10-100 ng genomic DNA, and 1 U of Hot FIREpol DNA polymerase (Solis Biodyne, Tartu, Estonia); the thermocycler profile consisted of one cycle of 10 min at 95 • C followed by 35-40 cycles of 95 • C for 15 s, 60 • C for 30 s, and 72 • C for 1 min; and a final extension step of 72 • C for 60 min concluded the profile. Reactions were analyzed using capillary electrophoresis with an ABI 3100 × l Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) using five-dye spectral calibration. Allele binning and genotype scoring were performed using GeneMapper 4.0 (Applied Biosystems). All samples were genotyped for all markers in the three panels in at least two independent multiplex reactions; uniplex reactions were also performed for some few recalcitrant markers and/or samples.

Study
Microsatellites were analyzed using MICRO-CHECKER 2.2.3 [35] to detect signs of stutter and null alleles when using the total material. In addition, to search for specific loci being candidates under directional or balancing selection, compared to neutral expectations, all microsatellites were analyzed in Arlequin 3.5.2.2 [36].
In Fstat 2.9.4 [37], we estimated heterozygosity (H O ), expected heterozygosity (H E ), allele counts, and F-statistics for the 10 geographical locations sampled. GENEPOP 4.0.10 [38] was used to estimate the deviations from Hardy-Weinberg equilibrium and genotypic linkage disequilibrium using an exact (probability) test for locus and locations [39]. Correction for multiple testing was performed using the false discovery rate (FDR) in the FDR software [40] with the "graphically sharpened method." The results from these initial tests of the 26 microsatellite loci showed that 49 pairwise tests of the population and loci exhibited deviations; however, after a correction using the FDR method (threshold < 0.0027), only 13 combinations exhibited deviations from Hardy-Weinberg equilibrium. There was no systematic combination of loci and population, and the deviations were possibly due to null alleles. Out of the 2760 combinations of loci and populations, only 158 combinations indicated significant linkage disequilibrium (LD) before FDR correction, and no pairwise combinations were in significant LD after an FDR correction (threshold < 0.000058). Thus, after this initial microsatellite quality assessment, a set of 24 interpreted neutral microsatellites (Supplementary Table S1) were used in all further population genetic analyses.
F ST values among localities was calculated in FSTAT 2.9.4 [37]. In order to visualize the genetic differentiation of red fox individuals in the different study areas, we used the package Factoextra [41] in R [42] to construct a neighbor-joining tree with 1000 bootstraps. Genetic relationships among individuals and localities were analyzed in STRUCTURE 2.2.3 [43], revealing the most likely numbers of Hardy-Weinberg equilibrium clusters analyzed using an admixture model and correlated allele frequencies with 100,000 burn-in steps and 300,000 Markov chain Monte Carlo (MCMC) iterations and 10 replicates for each of the runs for K populations = 1-10. Three sets of analyses were performed, the first with and the second without the LOCPRIOR option (i.e., using geographical information of sampling localities to help STRUCTURE when low differentiation existed), and the third analysis was done in a hierarchical manner by removing one divergent area and then rerunning with the remaining samples (revealing the cumulative number of K groups). Evaluation of the most likely number of K was based on methods in Evanno, et al. [44] and STRUCTURE-HARVESTER, version 0.6.1 [45]. Here, we particularly compared the patterns of lnP (data) and ∆K used for the inferences. The visual graphics from the best STRUCTURE run (most likely K) was added to the neighbor-joining tree (described above).
To assign individuals based on both genetic and spatial information, the clustering software Geneland [46] was used for the four localities. Geneland offers the opportunity to both identify genetic clustering and to infer what influence spatial clustering has on the population structure [47]. Geneland identifies clusters based on an attempt to maximize Hardy-Weinberg and linkage equilibrium and uses MCMC to estimate the number of clusters (K). In addition, Geneland incorporates the spatial location of individuals, which can help to elucidate the influence of geographic spread on, e.g., gene flow. MCMC iterations were set to 1,000,000, uncertainty on the coordinate was set to 500 m, and thinning was set to 100. The uncorrelated model was run 10 times where K was set as minimum 1 and maximum 10 to ensure the correct number of K [48]. When K was determined, the model was run an additional 25 times with a set K to be able to choose the run with the highest log posterior probability of population membership for further analyses. Pixels were calculated with 500,000 iterations and 300 was chosen for the number of pixels.
To estimate the gene flow between the four study areas, the population size parameter Θ (4N Eµ ), was estimated, where N E is the effective population size and µ is the mutation rate per site, as well as the level of historical gene flow (M: migration rate m per generation divided by the mutation rate µ per generation) between the red fox populations of Sweden and Norway. This process used MIGRATE 3.0.3 [49,50] on microsatellite repeats with a Brownian microsatellite mutation model, default settings, and a maximum likelihood (ML)-strategy with mutation rates treated as constant for all loci. MIGRATE estimates migration rates based on coalescent theory (i.e., migration back to the most recent common ancestor). In addition, we estimated contemporary gene flow (m: immigration rate per generation) between the Swedish and Norwegian red foxes using BAYESASS 1.3 [51] on microsatellite repeats. This software calculates migration rates over the last few generations based on the proportion of immigrants (m ij is the fraction of individuals in population i that are immigrants from population j). Since the mutation rates are unknown in the red fox, we used the default settings in MIGRATE and BAYESASS (two independent runs with 100,000 MCMC steps, where the first 1000 steps were burn-in, were considered sufficient to obtain convergence).

Results
The results from the genetic diversity analyses showed that the mean number of alleles per locus varied between 7.08 and 9.90, where Sweden North had the lowest number of alleles and Hedmark South had the highest number. Hedmark South also had the highest mean number of private alleles (0.71), whilst Sweden South had the lowest number of private alleles (0.34). Allelic richness varied from 6.67 to 5.73, with Hedmark South exhibiting the highest values and Sweden South exhibiting the lowest. Sweden South exhibited the lowest inbreeding coefficient (F IS ) with a value of 0.07, while the highest inbreeding coefficient was in Hedmark South with a value of 0.09 (Supplementary Table S2).
When testing the genetic differentiation among the four study areas, the pairwise genetic distance (F ST ) ranged from 0.008 to 0.208. Out of the six pairwise combinations, five pairwise comparisons were significantly differentiated. The two Swedish study areas were significantly differentiated from the two Norwegian study areas, as well as being significantly differentiated from each other. However, the two Norwegian study areas did not exhibit significant differences in F ST (Table 2). Table 2. Pairwise comparisons of genetic distance (F ST ) among the four study areas of red foxes. Significant differentiation after adjustment for multiple comparisons is marked in bold. The population genetic analysis revealed that the most likely number of clusters estimated by Structure and Structure Harvester was two (K = 2, Figure S2). Assigned clusters indicated that the first cluster consisted of the two Norwegian study areas, as well as the Sweden North study area. The Sweden South study area was assigned to a separate cluster. No further sub-structuring was found when removing the Sweden South study area and re-running Structure ( Figure 2).
The neighbor-joining tree exhibited a main divide with moderately well-supported bootstrap values (0.67) into two separate branches, where the first branch consisted of the two Swedish study areas, and the second branch consisted of the two Norwegian study areas ( Figure 2).
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 14 the highest inbreeding coefficient was in Hedmark South with a value of 0.09 (Supplementary Table  S2). When testing the genetic differentiation among the four study areas, the pairwise genetic distance (FST) ranged from 0.008 to 0.208. Out of the six pairwise combinations, five pairwise comparisons were significantly differentiated. The two Swedish study areas were significantly differentiated from the two Norwegian study areas, as well as being significantly differentiated from each other. However, the two Norwegian study areas did not exhibit significant differences in FST ( Table 2). Table 2. Pairwise comparisons of genetic distance (FST) among the four study areas of red foxes. Significant differentiation after adjustment for multiple comparisons is marked in bold.

0.008
The population genetic analysis revealed that the most likely number of clusters estimated by Structure and Structure Harvester was two (K = 2, Figure S2). Assigned clusters indicated that the first cluster consisted of the two Norwegian study areas, as well as the Sweden North study area. The Sweden South study area was assigned to a separate cluster. No further sub-structuring was found when removing the Sweden South study area and re-running Structure (Figure 2).
The neighbor-joining tree exhibited a main divide with moderately well-supported bootstrap values (0.67) into two separate branches, where the first branch consisted of the two Swedish study areas, and the second branch consisted of the two Norwegian study areas (Figure 2). Geneland reached a different conclusion regarding the number of population genetic clusters as it inferred the most likely number of clusters to be 1. The individual red foxes were spread across the geographical landscape, though they were still assigned to the same genetic cluster (Supplementary Figure S1). All individuals were assigned to a cluster, and there was no evidence of ghost clusters (i.e., clusters that contained no individuals).
The historical gene flow analysis resulted in almost identical Θ values for the two Swedish and the two Norwegian study areas. Θ values were calculated to be the following: Hedmark South and Hedmark North = 0.098, Sweden North and Sweden South = 0.097. In total, only one of the Geneland reached a different conclusion regarding the number of population genetic clusters as it inferred the most likely number of clusters to be 1. The individual red foxes were spread across the geographical landscape, though they were still assigned to the same genetic cluster (Supplementary Figure S1). All individuals were assigned to a cluster, and there was no evidence of ghost clusters (i.e., clusters that contained no individuals).
The historical gene flow analysis resulted in almost identical Θ values for the two Swedish and the two Norwegian study areas. Θ values were calculated to be the following: Hedmark South and Hedmark North = 0.098, Sweden North and Sweden South = 0.097. In total, only one of the comparisons of historic migration rates were significantly different based on 95% confidence intervals. The only significant difference in migration rate was found between Hedmark North and Sweden North, with a significantly higher migration from the Norwegian to the Swedish study area. In more detail, the rate was highest from Hedmark North to  (Table 3). Table 3. Coalescent estimates of historic mutation-scaled population sizes (Θ = 4N E µ) and migration rates (M ij = m ij /µ) with standard error between the four study areas from Norway and Sweden, estimated with MIGRATE (M). Significant migration rates are marked in bold. The inferred contemporary migration rates from BAYESASS (fraction of population i that originates from population j) indicated a significantly higher migration from Sweden North to Hedmark North with a rate of m = 0.30 (CI: 0.28-0.32), while the current migration from Hedmark North to Sweden North was estimated to be m = 0.007 (CI: 0.0003-0.0014, Figure 3, Table 4). There was also a significant degree of migration from Sweden South to Hedmark North with a value of m = 0.30 (CI: 0.28-0.32), and from Hedmark South to Hedmark North with m = 0.31 (CI: 0.30-0.32). All other migration rates were small and not significantly different, with a range from 0.006 (Hedmark South to Sweden South) to 0.01 (Sweden North to Hedmark South).

Discussion
Red foxes in Norway and Sweden exhibited limited differences in genetic diversity. The pairwise comparison of genetic differentiation (FST) indicated that all study areas exhibited low to moderate differences in genetic distance from each other, except for the two Norwegian study areas. Similar results were found by Norén et al. [7] in their study of the red fox in Fennoscandia, where low to moderate FST values were found, indicating high degrees of gene flow across Fennoscandia. The neighbor-joining tree separated the two Swedish study areas from each other and from the Norwegian study areas, but the two Norwegian study areas were grouped close to each other, though with limited bootstrap support. Analyses for the population genetic structure in the Norwegian and Swedish red foxes suggested limited genetic structuring, with high amounts of mixing between localities. The initial STRUCTURE analysis suggested that the four study areas consisted of two separate clusters, with Sweden South clustering separately from the remainder of the locations. Geneland, on the other hand, suggested only one cluster. Our results of limited genetic differentiation and population genetic structure is therefore in line with previous studies that indicated limited genetic structure in red fox populations (e.g., References [52,53]).
Several studies have investigated population genetic structure in wild red fox populations, though only a handful of studies revealed signs of a strong genetic structuring (e.g., References [54,55]). This is most likely due to the red fox's ability to traverse long distances, and thereby exhibit a high degree of gene flow [7]. However, Frati et al. [54] found that the red fox in the Mediterranean basin were of two distinct clusters. This distinct clustering was allocated to the presence of two putatively different colonization waves during quaternary glaciation. Cohen et al. [55] found evidence of four population genetic clusters in Israeli red foxes using microsatellite markers. Here, the genetic structuring was attributed to anthropogenic limits to connectivity. Likewise, Wandeler et al. [56], using microsatellite markers, observed weak population genetic differentiation in red foxes in the city of Zürich, Switzerland, which was attributed to a local river barrier.

Discussion
Red foxes in Norway and Sweden exhibited limited differences in genetic diversity. The pairwise comparison of genetic differentiation (F ST ) indicated that all study areas exhibited low to moderate differences in genetic distance from each other, except for the two Norwegian study areas. Similar results were found by Norén et al. [7] in their study of the red fox in Fennoscandia, where low to moderate F ST values were found, indicating high degrees of gene flow across Fennoscandia. The neighbor-joining tree separated the two Swedish study areas from each other and from the Norwegian study areas, but the two Norwegian study areas were grouped close to each other, though with limited bootstrap support. Analyses for the population genetic structure in the Norwegian and Swedish red foxes suggested limited genetic structuring, with high amounts of mixing between localities. The initial STRUCTURE analysis suggested that the four study areas consisted of two separate clusters, with Sweden South clustering separately from the remainder of the locations. Geneland, on the other hand, suggested only one cluster. Our results of limited genetic differentiation and population genetic structure is therefore in line with previous studies that indicated limited genetic structure in red fox populations (e.g., References [52,53]).
Several studies have investigated population genetic structure in wild red fox populations, though only a handful of studies revealed signs of a strong genetic structuring (e.g., References [54,55]). This is most likely due to the red fox's ability to traverse long distances, and thereby exhibit a high degree of gene flow [7]. However, Frati et al. [54] found that the red fox in the Mediterranean basin were of two distinct clusters. This distinct clustering was allocated to the presence of two putatively different colonization waves during quaternary glaciation. Cohen et al. [55] found evidence of four population genetic clusters in Israeli red foxes using microsatellite markers. Here, the genetic structuring was attributed to anthropogenic limits to connectivity. Likewise, Wandeler et al. [56], using microsatellite markers, observed weak population genetic differentiation in red foxes in the city of Zürich, Switzerland, which was attributed to a local river barrier.
Teacher, et al. [52], however, detected some signals of isolation by distance in their study of phylogeographic patterns using mtDNA markers, but they did not find any evidence of a phylogeographic structure of the red fox in Europe. Likewise, Gachot-Neveu et al. [57], in their study using random amplified polymorphic DNA (RAPD) markers, did not detect any genetic clustering in a rural red fox population in France, most likely due to high levels of gene flow. Population genetic structuring was not detected in six areas in northeastern Poland using microsatellite markers. Here, Mullins et al. [58] concluded that genetic differentiation among red fox populations in their study was low due to a high connectivity and gene flow.
Norén et al. [7] identified genetic differentiation among red foxes from southern, western, and northern areas of Fennoscandia, and the foxes exhibited a high connectivity and asymmetric migration rates. The authors [7] also revealed that the historical genetic differentiation was much stronger than the contemporary differentiation, which indicates that gene flow has increased in the Fennoscandian red fox population. Our results agree with the previous study in that red fox exhibited a high degree of genetic diversity and gene flow due to a large dispersal ability.
In summary, our results demonstrate that there was a high amount of gene flow in Scandinavian red foxes, probably due to frequent movements across long distances and across country borders, though with some limited differences in local populations. The extent of gene flow entails that Norwegian and Swedish populations are near identical in genetic structure and can be viewed as one, large, panmixing population.
A previous study showed a general trend of a northward dispersal patterns in Fennoscandian red foxes [7]. In our study, we revealed contrasting patterns between historic and contemporary gene flow. We observed a high degree of historic gene flow from Hedmark North to Sweden North, whereas contemporary gene flow followed an opposite pattern, i.e., from south to north. The highest gene flow was estimated to be from Hedmark South to Hedmark North, as well as from Sweden South and North to Hedmark North. Dispersal in a northern direction was expected from climate-induced range shifts (e.g., References [59,60]). However, the weak genetic structuring and high amount of gene flow between the geographical regions also indicated that gene flow occurred from several different source populations, and not only in a northern direction, which may indicate range expansion [7]. Changes in environmental conditions, or biotic changes, such as prey availability, may induce population increases, which again can influence the dispersal and expansion from local areas [61]. With increasing temperatures due to climate change, in addition to other anthropogenic changes, has increased the range of suitable red fox habitats northwards [62,63]. Another influential factor may be the major outbreak of the parasitic disease sarcoptic mange (Sarcoptes scabiei) that took place in 1975 with its origin in North-Central Sweden. The outbreak resulted in a population decrease of up to 90% in the affected red fox populations [64]. This reduction in red fox density may have created more vacant territories, inducing recent recolonization due to a population recovery [65]. Our findings are in line with the findings of Norén et al. [7] in their study of the geographic patterns of expansion in the Fennoscandian red fox. By comparing historical and contemporary samples of red foxes from Fennoscandia, they revealed low genetic structuring, high connectivity, and asymmetric migration rates, which indicates source-sink dynamics in the Fennoscandian red fox population. The same pattern was found in our study, and this indicates that the weak genetic structuring and high amount of gene flow is in line with the classic signs of a range expansion or recolonization. It is likely that the major outbreak of sarcoptic mange created sink areas with a low red fox population density, and that these sink areas are now being recolonized from a range of sources located in a southerly direction. It is possible that the shifting direction in migration rates from the historical southerly direction to a more northwesterly direction in the contemporary analyses in this study was due to historic large-scale biotic changes that has induced a range shift or changed the direction of the population expansion. As the outbreak of mange in 1975 started in North-Central Sweden [64], it is a possibility that the outbreak first created sink areas available for recolonization in a southern direction before the outbreak spread northwards. However, to investigate this, thoroughly historic samples of red foxes should be analyzed. Walton et al. [10] in their study of GPS-marked red foxes in Sweden and Norway found a statistically non-significant northwesterly trend in long-distance dispersal movements of red foxes that coincides with the findings in this study. Similar patterns of a preference for long-range dispersal in a northerly direction has been described for several other red fox populations [3,66,67]. Disregarding the potential effect of sarcoptic mange to create sink areas, northern areas may in themselves attract immigration due to lower population densities and thereby more available territories. However, in general, productivity decreases with latitude, which again may influence the tolerance level for dense red fox populations [68]. Regardless, the vast distances red foxes have been observed to disperse, in combination with the evidence of dispersal in a northwesterly direction, reveals a high range expansion in a northwesterly direction that can influence the spread of parasites and disease into Norway [3,10,63].

Considerations for Management Authorities
In our study, there was a high amount of gene flow between red foxes across the border between Norway and Sweden. Among our study sites, we have included red foxes from three areas known to have occurrences of E. multilocularis and our results indicate that a large degree of red fox dispersal took place from these localities and across the border to Norway. Echinococcus multilocularis is a highly dangerous zoonotic parasite that is considered to be an advancing disease in Europe and can be fatal for humans [69,70]. The spread of the parasite into new areas poses serious concerns for the public harvesting of, e.g., bilberry (Vaccinium myrtillus), and thus human health [71]. E. multilocularis has not yet been detected in Norway, but in 2010, the first infected red fox individual was detected on the west coast of Sweden [23]. Since then, positive tests for the parasite has been found in four additional areas in Sweden, and five areas are now known to be infected, namely Borlänge, Nyköping, Katrineholm, Vetlanda/Växjö, and Uddevalla [24].
The Norwegian authorities have implemented a mandatory vaccine program for the import of dogs and cats across the border from Sweden as a stage in controlling the spread of the parasite [72]; however, there is no program implemented to decrease the spread of the parasite through wildlife. The results from our study suggest that it may not take too long before E. multilocularis will occur in Norway, or it may already be present, but not yet detected.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/9/24/5289/s1. Figure S1: Spatial Bayesian clustering for the four study areas from the software Geneland. Individual foxes are represented as black circles. A) Spatial clustering into one cluster, B) posterior probability map of belonging to cluster 1. Figure S2: ∆K value for the most likely number of clusters from Structure Harvester. Table S1: Description of the designed primer sequences and multiplex panels, including primer concentrations and labels, observed allele numbers and size range, and the expected allele size deviations compared to using established primer pairs previously utilized in single-plex or multiplex reactions with fewer markers. Table S2