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.
Introduction
Genetic structures are generally the outcome of interactions among species' social traits, geographic features limiting dispersal, and historical distribution patterns (Wright 1943; Slarkin 1985; Hartl and Clark 1997; Epperson 2003), although ecological processes have also been shown to play an important role (Pilot et al. 2006). Information on the genetic structures of populations are therefore relevant for a broad range of topics ranging from evolutionary biology to directly applied environmental resource management (Epperson 2003; Manel et al. 2003). Genetic structure is particularly relevant considering the current rate of broad scale environmental change, largely driven by climate change and landscape transformation (Nelson 2005). Such environmental change can alter distribution patterns of species at a range of scales (Parmesan 2006; Sexton et al. 2009), and subsequently also the genetic structures of populations (Alsos et al. 2012; Franks and Hoffman 2012; Theodoridis et al. 2020). Furthermore, genetic variation is the raw material for evolution and therefore crucial for adapting to changing environmental conditions (Barrett and Schluter 2008).
Southern Africa may suffer significantly from projected global warming, primarily related to increased frequency and severity of droughts (Makholela et al. 2017). The direct ecological consequences of a potential increase in the magnitude and frequency of drought, induced by global warming, are difficult to predict in southern African arid grasslands and savannas (Le Houérou 1996; Sankaran 2019). However, increased woody plant cover, coupled with accelerated spatial expansion and intensification of agricultural activities are likely to have substantial ecological effects (Dobrovolski et al. 2011; Midgley and Bond 2015; however, see Hoffman et al. 2019 for an alternate view), including range shifts, range contractions and species loss of terrestrial mammals (Thuiller et al. 2006; Visconti et al. 2011). This expected range contraction and fragmentation is likely to influence interpopulational genetic exchange as well as genetic diversity within populations of southern African mammals.
The bat-eared fox (Otocyon megalotis) is a small, primarily insectivorous canid that has a disjunct distribution in southern and eastern Africa (Clark Jr 2005). The species is mainly confined to open grasslands and largely relies on termites (Isoptera) and other arthropods for its diet (Klare et al. 2011; Jumbam et al. 2019). Although classed as ‘Least Concerned’ in a recent South African IUCN Redlist assessment (Dalerum et al. 2016), quantification of spatial population structures and genetic diversity were identified 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).
Materials and methods
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 (Dalerum et al. 2017). 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/2009, 2536/2009, 2571/2010,2572/2010; Kalahari: 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 MgCl2 and 0.35 µl of 2.5 units of HotStarTaq polymerase (Qiagen, Hilden, Germany) and distilled H2O. 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).
Data analysis
Analysis of mitochondrial data
We used Arlequin (version 3.5, Excoffier et al. 2005) to calculate gene diversity, nucleotide diversity and population pairwise FST. 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 FIS values were calculated in Arlequin 3.5 (Excoffier et al. 2005); where significant, positive FIS values indicate a higher relatedness than expected from random mating.
Genetic differentiation was quantified as population pairwise FST 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).
Results
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 FST = –0.026 (p = 0.640).
Microsatellite data
Microsatellite multilocus genotypes were obtained for 24 bat-eared foxes. Microsatellite data are available from FigShare at http://doi.org/10.6084/m9.figshare.14769201. All six loci were polymorphic with observed heterozygosity levels ranging from 0.42 to 1.00 (Table 1). Significant heterozygote deficiencies were detected in two loci in Benfontein and three loci in the Kalahari (Table 1). There was no significant difference in allelic richness (χ2 = 0.026, df = 1, p = 0.871) or individual multilocus heterozygosity (χ2 = 0.001, df = 1, p = 0.976) between Benfontein and the Kalahari. Population-specific FIS values across all loci were positive with FIS = 0.137 (p = 0.011) in Benfontein and FIS = 0.254 (p < 0.001) in the Kalahari.
There was a significant divergence between sampling regions (FST = 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 (FIS = 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 (tafz_a_1942204_sm2806.pdf)). 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).
Figure 2:
Minimum spanning network for five haplotypes of the cytochrome b (406 bp) found in 19 bat-eared foxes (Otocyon megalotis) from Benfontein and the Kalahari, South Africa. The size of each respective circle corresponds to the total number of observed foxes with a specific haplotype, and the relative pie charts reflect the proportion of those foxes from each study area
Table 1:
Hardy–Weinberg proportions for six microsatellite loci in Benfontein and the Kalahari
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 (Dalerum et al. 2017; 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.
Figure 3:
(a) Biplot of the coordinates of the first two dimensions from a Principal Component Analyses based on inter-individual genetic distance of South African bat-eared foxes (Otocyon megalotis). (b) Barplot giving the probability of membership in one of two a priori defined populations in a genetic assignment analysis
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
Although we observed only moderate relatedness among animals within each population, the observed positive and significant FIS 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 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.
Acknowledgements
We are grateful to the De Beers Ecology Division and the Diamond Route for permission to work in the Benfontein Game Reserve and to the staff and managers at Benfontein for logistic support. We similarly thank the Kalahari Research Trust, the Kuruman River Reserve and the Kalahari Meerkat Project for logistical support and access to the Kuruman River Reserve. The research was funded by the National Research Foundation of South Africa (90491 awarded to Fredrik Dalerum and TTK1206041007 awarded to ALR), as well as Spanish Ministry of Economy and Competitiveness funding (RYC-13-14662 awarded to FD). Joana M da Silva provided useful comments to an earlier draft of the manuscript.
