Genetic Variation between and within Two Populations of Bat-Eared Foxes (Otocyon megalotis Desmarest, 1822) in South Africa

Information on genetic variation within and among populations is relevant for a broad range of topics in biology. We use a combination of mitochondrial and nuclear microsatellite markers to evaluate genetic variation within and between two populations of bat-eared foxes (Otocyon megalotis Desmarest, 1822) in South Africa. The bat-eared fox is a small canid occurring in southern and eastern Africa. The species is currently not threatened with extinction, but a lack of information on genetic diversity has been identified as a deficit for its future conservation. We observed low to moderate genetic differentiation between the two geographically separated populations, but neither mitochondrial nor nuclear microsatellite markers suggested that there have been dispersal barriers between them. Similar genetic diversity within both populations was contrasted by interpopulational differences in relatedness variation among males and females. A high genetic relatedness within both populations, indicated by mitochondrial data, is likely caused by a common historical origin or a combination of species-specific social organization and environmental dispersal constraints. We call for further research on the genetic divergence of bat-eared fox populations as well as on the genetic consequences of interactions between environmental characteristics and social organization in this species.

as research priorities for the species. These priorities were identified because of a lack of existing information, which reduces our ability to predict the future consequences of ongoing environmental change (Midgley and Bond 2015).
Here we use mitochondrial and nuclear genetic markers to evaluate the population genetic structure of bat-eared foxes from two study areas in contrasting environments in South Africa. Separated by approximately 410 km, one of the areas lies in the transition between Karoo scrubland and arid grasslands, and the other in an arid Kalahari environment. In this pilot study, we simultaneously present assessments of bat-eared fox genetic structures, both within and between populations, as an addition to the previous scant information on genetic variation in this species (Wright et al. 2010;Kamler et al. 2013;Westbury et al. 2017).

Sample collection
Samples were collected as part of long-term projects on the behavioural ecology of bat-eared foxes (le Roux et al. 2014;Jumbam et al. 2019). The studies were carried out on Benfontein Game Reserve (Benfontein) and the Kuruman River Reserve (Kalahari), which lies in central (Benfontein, 28°99′ S, 24°81′ E) and northern (Kalahari, 26°99′ S, 24°81′ E) South Africa and are separated by approximately 410 km (Figure 1). Both reserves are privately owned, but differ in size (Benfontein 11 000 ha, Kuruman River Reserve 33 000 ha). Benfontein lies in a transitional vegetation zone consisting of dry Karoo, grassland and Kalahari thornveld (Schulze and McGee 1978). The area has a semiarid climate, with a dry season comprising March to August and a wet season from September to February. Main habitat types include open grasslands, low-growth scrubland and acacia-dominated open savanna . The study area in the Kalahari primarily consists of Kalahari thornveld vegetation and sparsely grassed sand dunes (Mucina and Rutherford 2011). The area experiences two distinct seasons, a cold-dry season (May-September) and a hot-wet season (October-April) (Russel et al. 2002).
Fourteen foxes were sampled on Benfontein between May 2008 and March 2013, and 10 foxes in the Kalahari during February and April 2014. The foxes at Benfontein were all captured within 5.8 km of each other, with nine foxes captured within 1.5 km of one another. At two other locations, three and two foxes were captured, respectively. The foxes in the Kalahari were all captured within 5.5 km of each other. Three of the foxes in the Kalahari were captured at the same location, with two foxes each at two other locations. On Benfontein, eight foxes were males and six females. In the Kalahari, five were males, four females and one was without information on sex. Small pieces of skin were cut from the ear-tips during routine captures and stored in 95% EtOH at −20 °C until the analyses were done. Captures at both sites were carried out using cage traps and subsequent immobilization with ketamine hydrochloride (2.2-5.5 mg kg −1 ) and medetomidine hydrochloride (44-111 µg kg −1 ), which was reversed with atipamezole (0.2-0.5 mg kg −1 ). The captures were carried out under permits from the Northern Cape Department of Nature Conservation (Benfontein: 2535(Benfontein: /2009(Benfontein: , 2536(Benfontein: /2009(Benfontein: , 2571(Benfontein: /2010(Benfontein: ,2572(Benfontein: /2010Kalahari: 476/2/2013) as well as permissions from the Animal Use and Care and Research Committees of the University of Pretoria (ec031-07 and V047/10) and the ethics committee from the University of the Free State (11/2013).

Genetic analyses
Amplification and sequencing of mitochondrial data DNA extractions were performed in a separate laboratory, isolated from post-PCR reaction products. A 406 bp fragment of the mitochondrial cytochrome b gene was amplified using primers L14724 (Meyer and Wilson 1990) and H15149 (Kocher et al. 1989). Each PCR reaction had a total volume of 25 µl, containing 2 µl of DNA template, 2.9 µl 10× PCR reaction buffer, 1.5 µl of primer, 2.3 µl of 10 mM dNTP, 2.9 µl of 25 mM MgCl 2 and 0.35 µl of 2.5 units of HotStarTaq polymerase (Qiagen, Hilden, Germany) and distilled H 2 O. The amplification was carried out on an MJ Research PTC-100 Thermal Cycler using the following cycling conditions: 95 °C for 10 min, followed by 35-40 cycles 94 °C for 1 min, 50-53 °C for 1 min and 72 °C for 1 min, and ending with a final elongation step of 72 °C for 5 min. PCR products and negative controls were analysed through gel electrophoresis at 75 V for 30 min. Successfully amplified products were sent to the company Macrogen Inc. (Amsterdam, Netherlands) for sequencing. All sequences were aligned and visually monitored using Bioedit (Hall 1999).

Amplification and sequencing of nuclear microsatellites
All samples were genotyped at six microsatellite loci (CXX140, CXX250 [Ostrander et al. 1993], CPH3 [Fredholm and Winterø 1995], 758, 606 [Ostrander et al. 1995] and 671 [Ostrander et al. 1995]), originally developed for the dog genome. A protocol described for the arctic fox was used (Hasselgren et al. 2018), including PCR reaction mixtures with a total volume of 10 µl per samples, of which 2 µl were DNA extract, 1 µl of 10× PCR reaction buffer, 1 µl of 10 mM dNTP, 0.15 µl of 2.5 units of HotStarTaq polymerase (Qiagen, Hilden, Germany) and water to reach the final volume of 10 µl. The amplification was carried out on a Veriti Thermal Cycler (Applied Biosystems, Foster City, CA, USA) using touchdown cycling conditions: 95 °C for 15 min, followed six cycles of 94 °C for 30 s, 58-52 °C for 30 s (Δ-1 °C per cycle), and 72 °C for 1 min, followed by 32 cycles of 94 °C for 30 s, 52 °C for 30 s, and 72 °C for 1 min, and ending with a final elongation step of 72 °C for 10 min. The PCR products were size determined using a LIZ-500 size standard (Thermo Fisher Scientific, Waltham, MA, USA) on an ABI3730 capillary sequencer (Applied Biosystems, Foster City, CA, USA) at Macrogen Inc. (Amsterdam, The Netherlands).

Analysis of mitochondrial data
We used Arlequin (version 3.5, Excoffier et al. 2005) to calculate gene diversity, nucleotide diversity and population pairwise F ST . from the mitochondrial sequences. To test for significant genetic differentiation between sampling localities, 1 000 permutations were done. The relationships between haplotypes were visualised using a minimum spanning network in PopArt (Leigh and Bryant 2015).

Analysis of nuclear microsatellite data
We used an exact test consisting of 100 000 dememorisation steps implemented in Arlequin (version 3.5) to examine deviations from Hardy-Weinberg proportions for microsatellite DNA. Allelic richness was calculated in FSTAT2.9.3.2 (Goudet 1995). Individual multilocus heterozygosity (MLH) was calculated as the proportion of heterozygous loci in relation to the number of successfully genotyped loci. Differences in allelic richness and multilocus heterozygosity were compared between the two sampling areas using a Kruskal Wallis test implemented in the statistical environment R (version 4.1.0, http://r-project.org). As an approximation of population inbreeding levels, population-specific F IS values were calculated in Arlequin 3.5 (Excoffier et al. 2005); where significant, positive F IS values indicate a higher relatedness than expected from random mating.
Genetic differentiation was quantified as population pairwise F ST with significance testing using 1 000 permutations. A hierarchical analysis of molecular variance (AMOVA) was conducted, as well as a PCA based on inter-individual genetic distance in GenAlEx 6 (Peakall and Smouse 2006). Furthermore, a Bayesian MCMC approach in STRUCTURE was used to investigate the occurrence of a potential substructure between the sampling localities (Pritchard et al. 2009). To identify the most likely number of clusters, we set K = 1-4 with 10 replicates for each K and no prior population information. The number of clusters were identified using the Evanno approach (Evanno et al. 2005), implemented in STRUCTURE Harvester (Earl and von Holdt 2012). For the MCMC, 100 000 burn-in steps were set, followed by one million MCMC iterations. To assess individual assignment, we set K = 2 and used the LOCPRIOR model with sampling location. We calculated relatedness based on our microsatellite data for all pairs of animals according to Queller and Goodnight (1998). We used two-sample permutation tests to contrast relatedness between pairs of animals from the same populations to pairs from different populations, as well as to compare relatedness within pairs of animals from Benfontein to relatedness within pairs of animals from the Kalahari. We similarly used a permutation based two-way ANOVA to evaluate if differences in relatedness between males and females differed between the two areas, and two-sample permutation tests to contrast relatedness within males and females. Permutation tests were carried out using functions in the user contributed packages coin (version 1.4-0, Hothorn et al. 2008) and permuco (version 1.1.0, Frossard and Renaud 2019) for the R environment (version 4.1.0, https://r-project.org).

Mitochondrial data
We successfully sequenced 19 samples for the 406 bp cytochrome b fragment. In total, we identified five haplotypes of which four were present in Benfontein and three in the Kalahari. Mitochondrial sequences are available in GenBank (https://www.ncbi.nlm.nih.gov) with accession numbers MZ555761-MZ555765. Gene diversity was 0.643 ± 0.184 for Benfontein and 0.564 ± 0.134 for the Kalahari. Nucleotide diversity was also comparable between the two areas, with 0.0041 ± 0.003 in Benfontein and 0.002 ± 0.001 in the Kalahari.
The minimum spanning network shows that all identified haplotypes were closely related and that two out of the five haplotypes, both with central positions in the network, occurred in both areas ( Figure 2). In addition to the two shared haplotypes, Benfontein displayed two unique low-frequency haplotypes and the Kalahari one unique low-frequency haplotype. Genetic differentiation calculated from the mitochondrial divergence was non-significant with F ST = −0.026 (p = 0.640).
There was a significant divergence between sampling regions (F ST = 0.083, p = 0.001), with 8% of the total variance being partitioned between populations, 17% of the variance being partitioned among individuals and 74% within individuals (F IS = 0.187, p = 0.001). In agreement with this, there was low overlap in the PCA scores of the two sampling regions ( Figure  3a). Based on the Evanno approach, K = 2 was identified as the most likely number of clusters (Supplementary material Figure S1). For a solution using K = 2 clusters, population membership values were >0.9 for foxes from Benfontein and >0.8 for foxes from the Kalahari (Figure 3b).
Average relatedness was significantly lower between pairs of foxes from different populations (mean relatedness = −0.15, SD = 0.23) as compared to foxes from the same population (mean relatedness = 0.03, SD = 0.26, Z = 5.71, p < 0.001). While there were no differences between the two populations in average relatedness (Z = 0.76, p = 0.449), the variation in relatedness among males and females differed between the two areas (F = 7.081, p = 0.002, Figure 4). Males were significantly more related than females in the Kalahari (Z = 2.83, p = 0.001), whereas there were no differences between males and females in Benfontein (Z = 0.96, p = 0.339).

Discussion
The nuclear microsatellite markers suggested a moderate level of genetic differentiation between the two populations. We interpret these data as support for some isolation by distance (Wright 1943), but without any direct dispersal barriers between these two relatively distanced areas. Similar results have been found in the brown hyaena (Parahyaena brunnea) across southern Africa (Westbury et al. 2018). However, the mitochondrial data showed high levels of overlap, which indicate a shared historical origin. Such an interpretation agrees with previous findings on two other mesocarnivores, the black-backed jackal (Canis mesomelas) and the caracal (Caracal caracal), where mitochondrial markers revealed very limited genetic structuring across South Africa (Tensen et al. 2019). Although human activities may provide significant barriers to interpopulational movement of South African carnivores (Swanepoel et al. 2013), small species such as the bat-eared fox may more easily disperse across human altered landscapes. The area between our two study sites consists mainly of land used for low intensity cattle farming (Thompson 2019), which should not hinder the dispersal of bat-eared foxes. Our results may also be interpreted in accordance with previous assessments of genetic variation in the bat-eared fox across smaller spatial scales (Kamler et al. 2013), suggesting that the observed structuring may be caused by processes not directly related to distance, but rather social interactions and mating behaviour. We therefore require data from a broader range of bat-eared fox populations, including ones separated by more intensively transformed land, to be able to fully evaluate how this species will be affected by a potentially increased future habitat fragmentation.
Our results suggested similar levels of genetic variation and relatedness within both study populations, although we can not fully rule out that this lack of difference could have been caused by our relatively low sample sizes. However, relatedness among males and females appeared to have differed between Benfontein and the Kalahari. Intra-population genetic variation and structure are generally believed to be a consequence of dispersal tendencies and social organization (Parreira and Chikhi 2015). Both of these characteristics are related to mating patterns and resource utilization (Emlen and Oring 1977;Macdonald 1983;Sandell 1989). Genetic variation within bat-eared fox populations has subsequently been linked to both foraging  and mating habits as well as sex-biased dispersal (Wright et al. 2010;Kamler et al. 2013). Our two study areas lie in contrasting environments in terms of abundance and temporal fluctuation of food resources Jumbam et al. 2019). We interpret our results as indications of possible genetic consequences of such resource variation. However, despite our similar sample sizes of males and females, we can not fully rule out that the observed differences in sex-specific relatedness were not a sample artifact. Intraspecific variation in social organization, driven by geographic variation in resource use, has been previously suggested for various carnivores (Tannerfeldt and Angerbjörn 1998;Dalerum et al. 2006;Hulsman et al. 2010), although such resource-driven variation in social structure may not necessarily result in contrasting genetic variation (Holekamp et al. 2012). We argue that more studies are needed to fully disentangle the relationships between the distribution and abundance of critical resources, social structure, and genetic variation in this species.
Although we observed only moderate relatedness among animals within each population, the observed positive and significant F IS values in our microsatellite data suggests possible inbreeding or social structuring. Genetic consequences of social structuring was further implicated in the Kalahari by different genetic relatedness among males and females. Genetic structure within bat-eared fox  Figure 4: Average relatedness (mean ± SD) of bat-eared foxes (Otocyon megalotis) within one of two areas, Benfontein and the Kalahari, as well between these areas. Relatedness values were calculated using the Queller and Goodnight (1998) method for all pairs of animals, as well as for females and males separately populations appear to be dictated by social monogamy, low levels of promiscuity, and sex specific natal philopatry (Wright et al. 2010;Kamler et al. 2013). We suggest that the observed relatedness in our study could have been caused by an interaction of these species-characteristics and environmental constraints in dispersal, related either to the distribution of food resources (Macdonald et al. 1983) or predation risk (Palomares and Caro 1999). We did not have specific information of group identity in our animals, but group membership on Benfontein appeared to have been fluent (F Dalerum, pers. obs.). However, social units on both of our study sites consisted of social groups ranging from two to five adults (Le Roux et al. 2014; A Le Roux, pers. obs.), which contrasts the mostly monogamous pairs previously observed on a neighbouring study site to Benfontein (Kamler et al. 2013). We acknowledge that our study was made on a small number of individuals from only two sites. In addition, we only used a limited number of genetic markers. Despite these caveats, our study provides the first empirical data to show some genetic population structuring across South Africa. However, to fully resolve the genetic structure of bat-eared foxes within its southern range, further studies are needed utilizing data on more individuals from a broader set of sites and a more extensive set of genetic markers. We argue that research on the genetic differentiation of bat-eared fox populations separated by highly modified human landscapes would be particularly informative, as would studies on the genetic consequences of interactions between resource availability, predation risk, and species specific mating and social characteristics.
To conclude, the nuclear genetic markers suggested moderate genetic differentiation between two populations of South African bat-eared foxes that were geographically separated, but did not have any substantial dispersal barriers between them. Similar genetic diversity and relatedness in both populations were contrasted by differences in relatedness variation among males and females. Data on mitochondrial markers indicated a high level of genetic overlap within both populations, which we attribute to a shared historical origin of the two populations, as well as a combination of species-specific characteristics in mating patterns and social organization potentially coupled with environmental constraints in dispersal.