Influence of climate and landscape on genetic differentiation of aardvarks (Orycteropus afer)

As climate change accelerates, assessing how climate shapes gene flow and neutral and adaptive genetic differentiation on landscapes is increasingly important. Aardvarks (Orycteropus afer) are ecologically important in sub‐Saharan Africa but are sensitive to human pressures and increasing aridity. We used individual, population and landscape genetic approaches to infer the influence of landscape, climate and potential adaptive differences on gene flow.


| INTRODUC TI ON
In the face of rapid environmental change, describing and evaluating the causes of differentiation across the range of a species are important components of informed conservation and management.
Differentiation may arise from genetic drift after isolation of populations due to distance, landscape resistance to movement or barriers (Manel & Holderegger, 2013), at widely differing time scales depending on population size (Epps & Keyghobadi, 2015).Such isolation increases the loss of genetic diversity and limits the spread of beneficial alleles that arise locally (e.g., Creech et al., 2017).Local adaptation (i.e., parallel local adaptation, Kawecki & Ebert, 2004), can also be an important driver of differentiation for populations inhabiting heterogeneous landscapes, where different environments favour different genetic variants.Gene flow may disrupt local adaptation (Kawecki & Ebert, 2004), but when selection is strong, may be restricted by the lower fitness of migratory or hybrid individuals (Orsini et al., 2013).The reduction in gene flow and resulting differentiation has been referred to as isolation by adaptation (Orsini et al., 2013) or isolation by environment (Bradburd et al., 2013).
Climate differences are a major cause of adaptive differentiation (Keller et al., 2013), but climate-and recent changes in climate-may also influence landscape resistance to gene flow and genetic drift through effects on population size or habitat suitability (e.g., Castillo et al., 2016;Creech et al., 2020).Teasing those processes apart in genetic analyses is particularly relevant for inferring the consequences of rapid climate change.Some studies have used genomic approaches to identify particular variants associated with environmental gradients such as precipitation or temperature (e.g., Jordan et al., 2017).Others have used landscape genetic approaches to test for genetic differences associated with differences in climate, particularly across elevational gradients (e.g., Funk et al., 2016;Sexton et al., 2016).Such studies offer the potential to differentiate isolation by the environment from isolation by distance (Wright, 1943), resistance (McRae, 2006), or barriers, by testing multiple models simultaneously and interpreting differentiation associated with site differences rather than intervening landscapes as potential evidence of isolation by adaptation (Van Buskirk & van Rensburg, 2020).For many species, however, the scales or conditions at which isolation by distance, resistance, barriers or environment arise are unclear.While the importance of maintaining connectivity and gene flow among populations is generally well recognized (Rudnick et al., 2012), how landscapes shape those processes as well as adaptive variation often remains unevaluated.
Aardvarks (Orycteropus afer), the sole extant representative of the Order Tubulidentata, are a striking example of a species that is widespread, ecologically important and well-recognized, yet one for which population fragmentation and the causes and degree of genetic differentiation remain unassessed.Aardvarks occur across sub-Saharan Africa, are large (40-100 kg), nocturnal and generally solitary (Shoshani et al., 1988); their burrows create important habitat for many other species (Pike & Mitchell, 2013;Whittington-Jones et al., 2011).Up to 18 subspecies have been described from morphology, but the validity of those subspecies is unclear (Lehmann, 2009;Shoshani et al., 1988).Dispersal, population structure and phylogeography remain uninvestigated.While not considered at high risk given their wide distribution, local and regional distribution is not always well understood (Kingdon, 1997).Aardvarks are killed by humans for food or agricultural damage, are suspected to have been locally extirpated in some areas and have been rated as highly vulnerable in a regional conservation assessment (Freitag & Van Jaarsveld, 1997).Aardvarks apparently occur over wide precipitation and elevation gradients, raising the potential for adaptive differences, but aardvarks in arid ecosystems are considered vulnerable to drought and increasing aridity expected with climate change (Rey et al., 2017;Taylor et al., 2018;Weyer et al., 2020).Determining how landscape and climate shape genetic differentiation at local, regional and larger scales would enhance our understanding of conservation needs for this ancient and taxonomically distinct species (Kuntner et al., 2011).Examining these factors and the potential for climate-related adaptive differences could facilitate assessing how climate shifts, fragmentation or anthropogenic exploitation will affect populations.
Here, we evaluate the degree and potential causes of population differentiation of aardvarks by investigating population genetic structure at multiple spatial scales.After developing a suite of variable microsatellite loci, we used non-invasive methods to survey for aardvarks and collect faecal samples across South Africa, in two locations in Eswatini and at one location in Kenya.We estimated and described genetic structure within and between sampling areas using model-based and model-free approaches with population-and individual-based analyses.We conducted individual-based landscape genetic analyses in one large protected area (Kruger National Park, South Africa; hereafter, local scale) and across South Africa and Eswatini (hereafter, regional scale) to evaluate how landscape affected gene flow and genetic structure over local (~350 km) and regional (~1600 km) spatial extents.
At the local scale, we tested hypotheses about the influences of distance and landscape on gene flow using a locally developed habitat model (Epps et al., 2021).At the regional scale, we tested two sets of model predictors representing hypotheses regarding the effects of landscape on genetic differentiation.The first set, here- and the apparent absence from the driest regions of sub-Saharan Africa (Kingdon, 1997), we hypothesized that resistance to dispersal increased and thus gene flow declined with lower precipitation across the study area.Given that predicted habitat use by aardvarks was highest at intermediate levels of precipitation in Kruger National Park (Epps et al., 2021), we also considered an alternate hypothesis that regional resistance was lowest across areas with intermediate annual precipitation.Because steep, complex topography acts as a barrier for many terrestrial or semi-fossorial mammals (Badgley, 2010), we also hypothesized that terrain ruggedness increased resistance to dispersal.We further considered that areas of high temperatures might limit gene flow or that areas of either particularly high or low temperatures might limit gene flow.
The second set of hypotheses addressed potential adaptive differences using predictors described hereafter as 'difference models'.
We hypothesized that aardvarks in areas with divergent elevations, maximum temperature and average precipitation are adapted to different climatic regimes.Because gene flow may be less prevalent among populations of animals adapted for different climates or habitat types (Orsini et al., 2013), we predicted that genetic differentiation among individuals would increase with differences in elevation, maximum temperature and average precipitation at sampling locations.Finally, we used repeated detections of individuals to make inferences about space use, tested for closely related individuals to estimate dispersal distances, and investigated how relatedness and genetic differentiation among individuals varied with distance.

| Genetic marker identification and design
Aardvark-specific genetic markers have not been described and DNA samples are difficult to obtain given the cryptic nature of the species.The only previously published study of intraspecific genetic variation in aardvarks used mtDNA sequences amplified from universal primers and animals from European zoos (Pohlova et al., 2014).We sought variable genetic markers suitable for detailed investigations of population genetic structure by non-invasive sampling methods such as faecal DNA and that enabled reliable identification of individuals.For those purposes, microsatellites appeared most appropriate.
We extracted potential candidate microsatellite loci from the Orycteropus afer draft genome (Di Palma et al., 2012) using the Genome-wide Microsatellite Analyzing Toward Application (GMATA) software version 2.01 (Wang & Wang, 2016).We set parameters in GMATA to search for microsatellite motifs of 2-5 base pairs in length with a minimum of 10 repeats per motif.GMATA contains an internal dependency for Primer3 v2.2.3 (Untergasser et al., 2012), allowing the user to design primers for the amplification of identified microsatellites.For primer design, we set a target product length of 100-250 base pairs and considered 100 base pairs of each flanking region for primer binding.All other parameters in both GMATA and Primer3 were set to default values.This combination of parameters yielded 124,970 putative microsatellite markers, 92,390 of which had unique primer sequences within each flanking region and were therefore amplifiable.
We selected 84 microsatellites with tri-or tetranucleotide repeat sequences, each on a distinct scaffold within the genome assembly and with expected amplicon sizes between 110 and 220 base pairs.We initially tested each locus individually with unlabelled primers and a gradient of annealing temperatures (57-62°C); loci that produced clean amplicons were then dye-labelled and screened on an ABI 3730 DNA analyser for variability in samples from across the study area and to estimate product size more accurately before multiplexing.Reagent concentrations and thermocycling conditions for single-locus tests follow those described in Appendix S1 in Supporting Information.

| Study site selection
We chose to evaluate population genetic structure in South Africa and Eswatini, as the region has many accessible protected areas of varying sizes.We included one site from east Africa to provide a continental-scale comparison of genetic structure to contrast with our regional investigation in South Africa and Eswatini.Our study employed both population-and individual-level analyses.
We collected most of our samples from protected areas, but also from several locations where permission was granted by local landowners.All protected areas surveyed were reported by managing agencies or personnel as having aardvarks present or detected in the past.We sampled 7 protected areas in the Western Finally, we sampled at the Mpala Research Center in Kenya, providing continental-scale inter-population comparisons at distances from ~2600 to 4250 km (Figure 1).

| Collection of non-invasive samples
We collected faecal samples as sources of DNA.Aardvarks typically bury their faeces, often in characteristic scrapes, but sometimes in loose soil where they are foraging or burrowing (see Epps et al., 2021 for a guide to sign identification).We located samples by searching protected areas or other locations by vehicle or on foot.We attempted to survey systematically along travel routes in accessible portions of each protected area, except for Kgalagadi Transfrontier Park, where we searched only the southern portion of the South African section.When searching by vehicle, we drove slowly (<20 km/h) while watching for characteristic signs of aardvark digging; if seen, we searched the area on foot.
We followed track lines of individual animals for up to 3 h, examining areas of loose dirt and scrapes or searched for fresh digging when tracking was impractical.Fresher samples were prioritized when multiple faecal deposits were found within ~250 m, to limit duplicate samples from the same individuals in the same locations, but we often collected multiple samples because genotyping failure rates were high (see Section 3).Samples were stored in paper bags and dried within 24 h using solar heat (vehicle windshield) or in pots of sand over campfires or stoves.After field collection was complete, samples were subjected to 75°C for 45 min to satisfy import permit requirements and then stored at room temperature.Samples were collected and exported under per-

| Genotyping, identification of individuals, Hardy-Weinberg equilibrium, linkage disequilibrium, genetic diversity and related individuals
We extracted DNA from faeces and genotyped samples at 19 microsatellite loci in 4 pre-PCR panels (Table 1), with a minimum of 3 PCR replicates for any samples retained for analysis (Appendix S1).After generating genotypes, before evaluating genetic marker behaviour, we inferred whether any samples likely originated from the same individual (Appendix S1).After identifying sets of duplicate samples, we removed all but the most complete genotype for all subsequent analyses except those for space use (see below).To test for Hardy-Weinberg proportions (HWP) and linkage disequilibrium (LD), we grouped samples into putative populations based on location (hereafter, sample populations).
Each protected area where we detected aardvarks was treated as a sample population, as were clusters of incidentally collected samples provided by private landowners if separated by <75 km.
We used GENEPOP v. 4.1 (Rousset, 2008) to conduct global tests for HWP by locus and by population for both heterozygote excess and deficiency, using a Bonferroni correction for the number of loci or populations, accepted default settings for numbers of iterations and tested for LD between each locus pair across all populations.
F I G U R E 1 Protected areas searched for aardvark faecal samples in South Africa, Eswatini and a single site in East Africa (insert, top left); locations of faecal samples collected in South Africa and Eswatini (main map, italicized names for population locations, blue pentagons for individual sample locations).Black polygons (main map) indicate individual or aggregations of samples outside of protected areas used for population-level analyses, including Roggeveld, Victoria West, Prieska and Upington.

| Relatedness, population genetic structure and genetic assignment tests
We evaluated population structure in a multi-step process.First, we identified potential pairs or sets of close relatives (parent-offspring pairs or full siblings) and removed all but one of each set.
Estimating relatedness requires an accurate assessment of allele frequencies in each population, which requires avoiding including samples from strongly divergent gene pools while maximizing sample size.To that end, we conducted a preliminary analysis using Program STRUCTURE (Pritchard et al., 2000) to examine which sample populations could be combined.We conducted a 50,000-iteration burn-in followed by 100,000 sampling iterations for 5 replicate runs varying the number of clusters (K) from 1 to 5 as well as the use of informed priors and initiation at population identity.
We evaluated likely cluster number using the high assignment criterion (Epps et al., 2018), where the highest K for which every cluster had at least one individual assigned at some high probability (q) is selected; here, because of some low sample sizes, we employed q > 0.80.Subsequent groupings of samples are referred to hereafter as 'combined populations'.
After combining populations given preliminary population assignments, we used COLONY (Jones & Wang, 2010) to test for close relatives within combined populations, using the 18 loci exhibiting HWP proportions (see Section 3).We tested for full siblings using a full likelihood approach and a medium run length with a 'weak' prior for sibship, while allowing for polygamous males and females, considering all genotypes as 'offspring', and not scaling for sibship.For any full sibling or parent/offspring pairs or triads identified at probability >0.80 in COLONY or by visual inspection, we removed all but one individual from the dataset prior to further analysis of genetic structure.Then, we conducted final population assignment tests using Program STRUCTURE, using 200,000 iterations for burn-in, 800,000 iterations for sampling, and 5 replicate runs for each value of K from 1 to 5, using informed priors and initiation at population identity based on combined population membership, but excluding population identity for samples from under-sampled gene pools (Upington, Prieska, Eswatini).
Next, we evaluated genetic structure using both model-based and model-free approaches at population and individual levels.
After restricting the data set to the 18 loci that demonstrated HWP and removing all but one of any set of closely related individuals (see Section 3), we used GENEPOP to estimate inbreeding coefficients (F IS ), observed (H o ) and expected heterozygosity (H e ), population pairwise F ST values among sample populations and among combined populations; we used Arlequin v. 3.11 (Schneider et al., 2000) to test for significance of F ST values.We used bootstrapping over loci to estimate 95% confidence intervals for F ST and F IS values, performed using the hierfstat package version 0.5-11 (Goudet et al., 2022).
TA B L E 1 Primer sequences, repeat motif and number of repetitions in the reference genome (Di Palma et al., 2012) for microsatellite loci of aardvark, Orycteropus afer, and dye labelling information for four multiplex PCR panels, with allele size ranges for samples collected in South Africa, Eswatini and Kenya.We made inferences about the distances at which dispersal occurred by two methods.First, we measured the distances between closely related individuals identified using COLONY.Secondly, we made inferences about genetic neighbourhood size using correlograms (Borcard & Legendre, 2012).We estimated relatedness among individuals using Lynch and Ritland's (1999) estimator, calculated using SPAGEDI (Hardy & Vekemans, 2002), but restricted estimates to within combined populations of >20 samples.We generated a correlogram of relatedness estimates against the distance between individuals within each combined population, evaluating correlation at 10 distance intervals balanced to provide the same number of pairs in each category, to determine the scale at which relatedness stopped decreasing with distance.Finally, to infer how genetic neighbourhood size (Wright, 1946) varied with distance, we used the PCA-based estimates of genetic distance among all individuals (details in next section) within South Africa and Eswatini and evaluated correlation at 5 km intervals from 0 to 100 km, with greater distances collapsed into a single interval, to determine the scales at which individuals were more or less genetically similar than expected across the dataset.

| Individual-based landscape genetic analyses
We conducted individual-based landscape genetic analyses at two extents: regional (South Africa and Eswatini) and local (Kruger National Park).At the regional scale, we analysed how in Julia (version 1.6.3,Bezanson et al., 2012).We used the resistance predictors to address our hypotheses that steep, rugged We used linear mixed models (LMM) to assess the relationships between pairwise individual genetic distances and resistance or difference variables while controlling for individual identity in each pairwise comparison (Shirk et al., 2018).For genetic distance, we used pairwise Euclidian distances in a reduced dimension PCA space (Castillo et al., 2014;Shirk et al., 2010) where we used all 19 loci and retained enough axes to account for 95% of the variation in both regional (25 axes) and local scale (20 axes) analyses.To preserve local variation, instead of setting missing values to the centre of the principal component axes, we used the local mode (most common allele) for each locus to fill in missing data, drawing the mode from the appropriate combined population.Samples not falling into one of those (Eswatini, Prieska, Upington) were assigned a mode based on the global South Africa and Eswatini datasets.We also reanalysed using allele means (the typical approach, e.g., Jombart, 2008) to determine sensitivity to this choice.
To evaluate predictors of genetic distance, we used Akaike's Information Criterion (AIC) to select among univariate models, alternative variable transformations and a limited suite of multivariate models that included a difference and a resistance predictor (Burnham & Anderson, 2002).Predictors included direct or transformed applications of vector ruggedness (VRM), topographic position index (TPI), elevation, mean annual precipitation (MAP) and maximum temperature (T max ; see detailed explanations of those variables and their use as difference or resistance predictors in Table 2).We selected the model with the lowest AIC and highest AIC weight as the most supported models, but we considered all models within 2 AIC units of the top-ranking model competitive.We fit linear mixed models using the lme4 package (version 1.1-27.1;Bates et al., 2015) in R (version 4.1.1;R Core Team, 2021).We examined the Pearson correlation coefficients (r) among predictors using the basic R stats package, and avoided multivariate model structures with combinations of predictors with r > 0.7 (Table S1).We considered different transformations for some predictor values (Table 2), such as linear and quadratic transformations to describe resistance increasing with average TA B L E 2 Hypotheses and variables considered in regional (South Africa and Eswatini) and local-scale (Kruger National Park, South Africa) landscape genetic analyses.A DEM-derived measure of terrain ruggedness, 250 m resolution, defined as the difference between a pixel and its surrounding cells (Guisan et al., 1999;Jones et al., 2000;Weiss, 2001).d A DEM (Danielson & Gesch, 2011), 250 m resolution.Park.We chose this location as it had the largest spatial extent and genetic sample size (Table 3) and two models of predicted aardvark habitat suitability developed from non-invasive surveys across the park (Epps et al., 2021).For that analysis, we tested whether resistance surfaces created by inverting predicted suitability of the two habitat models (one fit using generalized linear models with a binomial error distribution and a logit link function, hereafter GLM, and the other fit using the maximum entropy estimator, hereafter ME, see Epps et al., 2021 for full details of developing those models), geographic distance or a null model best predicted genetic differentiation among individuals (Table 2).The habitat suitability models developed by Epps et al. (2021) estimated the relative use by aardvarks based on elevation, annual rainfall departure (i.e., difference from 550 mm, the optimum level of rainfall for termites in Kruger National Park), the normalized differential vegetation index, the percentage of sand in soil and water occurrence.

| Collection of non-invasive samples
From 2016 to 2018, we collected 253 aardvark faecal samples from 9 of the 11 protected areas in which we searched systematically (Figure 1; Table 3), as well as three focal areas and one additional sample (Prieska) outside of protected areas.We were unable to locate definitive signs of recent aardvark occupancy in Agulhas National Park or Kgalagadi Transfrontier Park, where we searched only the southernmost half of the South African section.

| Genetic marker identification, DNA extraction and genotyping
We developed primers for 19 variable microsatellite loci (Table 1).
Genotyping success was relatively low, presumably because most faecal samples were fully buried and mould was frequently observed on samples that appeared >24 h old, particularly in wetter environments.We attempted to genotype 196 samples.Of those, 125 yielded genotypes for at least 1 locus, with a mean proportion of 0.79 loci genotyped.

| Identification of individuals, Hardy-Weinberg equilibrium, linkage disequilibrium and genetic diversity
We determined that the 8 least informative loci met our standards for individual identification, yielding maximum P ID estimates of 1.82 × 10 −5 and P IDsibs of 0.0081 (Table S2).Therefore, we limited subsequent analyses to the 104 samples successfully genotyped at ≥8 loci, of which 47 yielded complete 19-locus genotypes.We identified 9 pairs of samples likely to originate from the same individuals and 5 instances where 3 samples appeared to originate from the same individual.Across all individuals, H e ranged from 0.118 to 0.863 (mean = 0.649) but exceeded 0.43 for all but one locus (Table S2).
Samples were collected from 12 primary areas: 9 protected areas and 3 clusters of samples in focal areas outside of protected areas separated by <64 km.Genotyping success from two sites in Eswatini was very poor, yielding only 1 sample with data for ≥8 loci and several other sites also had low sample sizes (Table 3).Therefore, tests for HWP and LD were conducted on the 10 sample populations with at least 2 samples.No evidence of heterozygote excess was observed; but heterozygote deficiency was detected in locus MK77186 (p < .0001,p critical = .0026).For sample populations, only Kruger National Park exhibited heterozygote deficiency (p = .0002;p critical = .005),likely because of genetic substructure caused by isolation by distance.We found no evidence of consistent LD among loci.Except for model-free approaches (e.g., PCA), subsequent analyses excluded MK77186 and thus were limited to the 18 loci not showing departure from HWP.

| Related individuals, population genetic structure and genetic assignment tests
Individual-based assignment tests using STRUCTURE (18 loci) suggested that 4 clusters were diagnosable (Figures S1-S3).Aardvarks in TA B L E 3 Areas where non-invasive genetic sampling for aardvarks was conducted, with numbers of samples collected, attempted to genotype and genotyped at ≥8 loci.Circles are linked to actual genotype location by grey lines where overlap would otherwise occur.Mean likelihoods of the data from replicated runs at K = 1-5 for this analysis are depicted in Figure S3.
Kenya and Kruger were strongly distinguishable, as were aardvarks in central South Africa (Tankwa Karoo, Roggeveld, Karoo, Camdeboo and Victoria West; Figures 2 and S2).Aardvarks in western/northern South Africa (Namaqua, Augrabies) formed a somewhat less distinct cluster, and samples from Upington, Prieska and Eswatini showed evidence of membership to gene pools not well represented by our sampling (Figures 2 and S2).
The four clusters identified using STRUCTURE were treated as combined populations for analysis and removal of closely related individuals (Table S3).Using COLONY and visual inspection of genotypes in areas with sample sizes <10, we identified two likely pairs in Kruger, no such pairs in Augrabies-Namaqua and in the central South Africa cluster, 3 closely related individuals in Tankwa Karoo and one pair in the Roggeveld area (Table S3).
After the removal of all but one of each closely related group of individuals, estimates of genetic diversity and genetic structure among samples clustered by protected and focal areas, as well as combined populations, showed considerable variation in genetic diversity.We also observed moderate to high genetic differentiation across South Africa and between eastern and South Africa (Tables 4 and S6), although small sample sizes for some sample populations added uncertainty.Moderate inbreeding coefficients (point estimates F IS from 0.05 to 0.15, 95% CIs did not overlap 0) were observed in Tankwa Karoo, Kruger and the Victoria West focal area, suggesting additional genetic structure or additional undiagnosed closely related individuals existed within those sampling clusters (Table S6).

| Inferences on space use and dispersal
Genotype matching demonstrated that aardvarks moved up to 7.3 km over the period during which faecal samples yielding usable DNA persisted (likely <1 month), suggesting maximum potential home ranges of 41 km 2 (circular) or 21 km 2 (elliptical).The median detection distance was only 1.6 km (Table S4).Full siblings or parent/offspring pairs detected were separated by 0.1-44.1 km (median = 8.2 km; Table S3).

| Individual-based landscape genetic analyses
Regional analysis revealed that both resistance and difference measures predicted individual genetic differentiation across TA B L E 4 Genetic diversity (observed [H o ] and expected [H e ] heterozygosity), inbreeding coefficients and genetic structure (F ST ) estimates using 18 microsatellite loci for combined populations of aardvarks, where sampling areas were combined on the basic of genetic similarity described by population assignment tests and after removal of 6 individuals determined to be closely related (full siblings/parent-offspring level) to individuals retained in the data set.between sampling sites (ΔT max ) strongly outperformed distance and all other resistance or differencebased measures in predicting individual genetic differentiation (Table 5).The distance-subsetting exercises for univariate models demonstrated that the inverse of mean annual precipitation resistance model performed well for datasets including most of the observations, including comparisons at close and far distances (Figure 3).Topographic ruggedness (VRM, Table 2) also appeared TA B L E 5 Landscape genetic results from individual-based regional scale analyses of aardvarks in South Africa and Eswatini.c Difference in average maximum temperature between sampling sites.

Population
d Difference in elevation between sampling sites.
e Vector ruggedness measure, where resistance increases with higher ruggedness.
f Difference in mean annual precipitation between sampling sites.
g Topographic position index, another measure of terrain ruggedness.
h Geographic distance between sampling sites.
i Gaussian transformation of mean annual precipitation, where resistance is highest at when mean annual precipitation is most different from the average.j Resistance increases with average monthly maximum temperature, 1970-2000.
k Gaussian transformation of elevation, where resistance is highest at relatively low and high elevations.
l Resistance increases as the square of the average daily maximum temperature.
m Gaussian transformation of average daily maximum temperature, where resistance is highest at low and high average maximum temperatures.
to perform well for datasets including all but the closest samples, although models based on lower resistance at moderate elevations performed well when analyses were restricted to the largest scales (samples >500 km apart).ΔT max operated most strongly at scales of <90 km but was not particularly effective as a univariate model (Figure 3; Table 5), indicating that its effects were most apparent when variation caused by aridity and distance in the intervening landscape was controlled for.Below upper thresholds of 80-500 km between samples, distance appeared to predict individual genetic differentiation more strongly than other models (Figure 3b), but models exploring lower thresholds did not support distance at any scale, suggesting that effects of distance were influenced by the closest pairwise comparisons (Figure 3a).Finally, the local scale analysis within Kruger National Park indicated that distance alone best predicted individual genetic differentiation on that landscape, compared to resistance measures from recently developed habitat models (Table S5), which aligns well with the scales at which distance was supported in the full dataset.In our individual-based landscape genetic models, our predictor for aridity, where resistance increased linearly as the inverse of mean annual precipitation, was present in the top three models (ΔAIC <5, Table 5).Thus, our hypothesis that arid areas act as a landscape barrier to aardvark gene flow was most strongly supported (Table 2).Lower rainfall likely results in less vegetation and thus fewer termites and ants on which aardvarks feed (Weyer et al., 2020)  another approach to addressing that important study component (Anderson et al., 2010).Scales of inference are often an artefact of study design, but because our study included close, separated and widely separated sampling areas, we were able to explore inferential resilience to variation in scale.For instance, exploring different upper thresholds demonstrated that our findings regarding the importance of aridity were stable at >500 km (Figure 3b), but a study focused on smaller scales might have missed this important relationship and the implications for climate change impacts on this keystone species.

| DISCUSS ION
Our individual-based analyses revealed the scale at which aardvark dispersal occurs.Across South Africa and Eswatini, a positive correlation of genotypes (Figure 4a) indicated dispersal typically occurs at <55 km.These estimates were supported by our detection of closely related individuals (likely full siblings or parent-offspring pairs) at distances up to 44 km (Table S3).Despite observing the longest distance between closely related pairs at Kruger National Park, average relatedness values there suggested that dispersal there may be most common over shorter distances (<21 km), whereas relatedness in central South Africa was higher than average at distances <56 km (Figure 4b).Thus, dispersal varies by landscape.Perhaps unsurprisingly (Keeley et al., 2017), our landscape genetic analyses within Kruger did not support resistance models, instead showing that distance best predicted individual genetic differentiation (Table S5).Kruger is a relatively homogenous landscape with no apparent barriers to aardvark dispersal, but perhaps territorial behaviours or high densities of large predators make longer-distance dispersal less common there.In other areas within our study, genetic discontinuities were suggested over small distances.For instance, samples near Upington reflected a different and largely unsampled gene pool compared to samples from Augrabies National Park <100 km distant (Figure 3).Likewise, in population-based analysis, large F IS values in Tankwa Karoo suggested substructure may exist within that arid park at the scale of <25 km (Table S6).  4 and S6, population based).
Population genetic structure between central and coastal South Africa also was weak, compared with much stronger differentiation between those locations and Kruger National Park in eastern South Africa, as well as Kenya (Table 4).As aardvarks are not limited to protected areas, this weak genetic structure in central South Africa suggests some continuity of aardvarks and stepwise gene flow across that region.We detected aardvarks in all the protected areas surveyed except Agulhas National Park and Kgaligadi Transfrontier Park (Figure 1), although South African National Parks (SANPARKS) listed aardvark in both areas.Kgalagadi is highly arid and both parks have sand as a dominate substrate, which may be less suitable for aardvarks.This highlights the need for more consistent surveys of aardvark presence or density even in protected areas.
Our microsatellite markers and methods for faecal DNA extraction offer new research opportunities.Non-invasive, systematic collection of faecal samples over time would allow detailed estimates of space use, relatedness, dispersal, population size or density and individual survival (e.g., Pfeiler et al., 2020Pfeiler et al., , 2021)), which have hitherto been very difficult to assess for aardvarks.For instance, our redetections of individual aardvarks suggested the potential for elliptical home range sizes of up to 21 km 2 (Table S4), in comparison to telemetry studies reporting home ranges of ~1.3-3.5 km 2 (Taylor & Skinner, 2003;Van Aarde et al., 1992).Thus, our study suggests that home range size in arid areas such as Tankwa Karoo National Park (Table S4), where prey densities likely are low, could well exceed reported values.Our markers showed moderate to high variability (Tables 4 and S4).However, we observed relatively poor amplification success.Amplification success for faecal DNA studies of other mammals can range from <25% in rainy or humid environments (e.g., Brinkman et al., 2010) to >95% in dry environments (e.g., Pfeiler et al., 2020).Future studies should prioritize collection of fresh or dry, exposed samples, explore other methods for rapidly drying or preserving samples or simply plan for low amplification success.
Our findings support recent research demonstrating that drought and temperature extremes negatively affect aardvark survival and thus aardvarks may be vulnerable to rapid climate warming (Rey et al., 2017;Taylor et al., 2018;Weyer et al., 2020).Aridity reduced gene flow, suggesting both habitat and dispersal limitations.
after 'resistance models', addressed hypotheses about how landscape between sampled areas might influence dispersal and thus gene flow.Given the vulnerability of aardvarks to drought and climate change (Rey et al., 2017; Taylor et al., 2018; Weyer et al., 2020) K E Y W O R D S climate, isolation by adaptation, landscape genetics, microsatellites, Orycteropus afer, Tubulidentata | 3 of 16 EPPS et al.

and
Northern Cape regions of South Africa (Agulhas National Park, Namaqua National Park, Augrabies National Park, Tankwa Karoo National Park, Karoo National Park, Camdeboo National Park), as well as the Kgalagadi Transfrontier Park on the border of South Africa, Botswana and Namibia and private lands in 4 areas in the Western and Northern Cape (Figure 1).Those areas were intended to provide inter-population comparisons at the ~75-1000 km scale.On the eastern side of southernmost Africa, we sampled Mlawula Game Reserve and Malolotja Nature Reserve in Eswatini and Kruger National Park, one of the largest protected areas in Africa.These sites provided the opportunity for comprehensive sampling and individual-scale analyses at distances of up to 360 km, inter-population comparisons within that region and inter-population comparisons across South Africa and Eswatini at the ~1000-1500 km scale.

2. 6 |
Individual-based inferences on space use and dispersalWe inferred space use by estimating geographic distances between locations where the same individuals were detected.For each individual, we used the longest distances between relocations (D) to create two indices of space use, including (1) a circular home range with a diameter of D, where area was estimated as * (D∕2) 2 and (2) an elliptical home range with a span D that represented the long axis and a span D/2 that represented the short axis, for which area was estimated as ( D ∕ 2) * (D ∕ 4).Because sampling often was constrained by proximity to roads or trails and redetections were sparse, we considered these estimates as indicators of the potential scale of home ranges, with the elliptical calculation representing a more conservative calculation.
terrain, higher elevations and areas of high aridity or low rainfall would constrain aardvark gene flow by impeding dispersal and limiting available habitat on intervening landscapes.Difference (Δ) predictors addressed the magnitude of difference between climate-associated measures including elevation, mean annual precipitation and average monthly maximum temperature assessed at individual genetic sample locations.Those predictors addressed our hypotheses that aardvarks living in habitats with widely differing climates would exhibit isolation by adaptation.We included geographic distance (i.e., the geodesic distance between two points on an ellipsoid, using R version 4.1.1,R Core Team, 2021 and the geosphere package version 1.5-14, Hijmans, 2021) and a null model as null hypotheses.
lower between aardvarks adapted to different climates Difference Elevation (ΔElev) Gene flow is lower between aardvarks adapted to different climates Difference Mean annual precipitation (ΔMAP) Gene flow is lower between aardvarks adapted to different climates Difference Null models were also considered.We described the hypothesized relationships between continuous landscape variables and resistance using four transformations that map continuous values into resistance values between r min = 1 and r max = 100.The four transformations are in the table footnotes where x is the value of the landscape variable, r is resistance, r max is the maximum resistance value, α is an exponential scaling factor, μ opt is the user defined value to centre the mean of the Gaussian transform and σ is the standard deviation of the Gaussian transform.a A digital-elevation-model (DEM)-derived measure of terrain ruggedness (Sappington et al., 2007), 250 m resolution, defined by a decomposition of slope and aspect into 3-dimensional vectors describing the terrain orientation of a grid cell neighbourhood.b We described a positive relationship between resistance values and landscape values using: r = indicating a linear relationship, for all transformations except for the squared T max transformation where α = 2, indicating a squared relationship.c

e
We described an expected relationship with lowest resistance at intermediate landscape values and increasing resistance at low and high resistance values using, r = (2017);Hijmans et al. (2005Hijmans et al. ( ), 1970Hijmans et al. ( -2000, 1 km resolution.gWe described a negative relationship between resistance and landscape values using: for all negative transformations, indicating linear relationships.h Fick & Hijmans (2017); Hijmans et al. (2005), 1970-2000, 1 km resolution.i Habitat suitability estimated by Epps et al. (2021) using two modelling approaches: Maxent and general linear modelling.maximum temperature, and linear and Gaussian transformations of resistance as a function of elevation, mean annual precipitation and average monthly maximum temperature.Gaussian transformations provide the lowest resistance at intermediate values (e.g., Castillo et al., 2014).Our circuit theoretic measures of resistance intrinsically included distance effects but we still used distance alone as a comparison model.As our sampling was by necessity clustered, we used two subsampling exercises (increasing and decreasing) to explore the dependency of our results on sample distribution.First, we pruned pairs of samples based on increasing intervening distances.We used 10 km intervals from 0 to 100 km, starting by removing all pairs separated by <10 km, then by <20 km.We pruned by 50 km intervals between 100 and 200 km, then by 100 km intervals out to the maximum extent of the dataset, assessing and plotting the AIC weight of each predictor variable fit in univariate models of genetic distance at each step while controlling for individual identity in each pairwise comparison.Second, we pruned sample pairs for decreasing distances by removing pairs in the largest distance intervals as above, with successive pruning down to the lowest distance intervals.These exercises allowed us to explore which predictors best explained genetic difference at scales ranging from individual dispersal distances to regional.At the local scale, we evaluated how landscape resistance shaped individual genetic differentiation within Kruger National The correlogram based on PCA distance among individuals for South Africa and Eswatini revealed significant (p < .05)genetic similarity among individuals separated by distances up to 55 km, although not at 30-40 km (Figure4a).The combined populations for central South Africa and Kruger National Park each contained enough individuals to justify the estimation of relatedness.Correlograms for each area showed that individuals were more related than expected out to 56 km in central South Africa and out to 21 km in Kruger National Park (Figure4b).
Samples from Upington (n = 2), Prieska (n = 1) and Eswatini (n = 1) were not included in this analysis due to the small sample size.95% confidence intervals estimated by bootstrapping over loci using the hierfstat package version 0.5-11(Goudet et al., 2022) are reported for F IS and population-pairwise F ST values.F ST values are reported in the 4 rightmost columns below the diagonal, with tests for significance performed in Arlequin v. 3.11 (Schneider et al., 2000) above the diagonal.a Includes Camdeboo NP, Karoo NP, Roggeveld, Tankwa Karoo NP and Victoria West.| 11of 16 EPPS et al.South Africa and Eswatini.Comparing analyses that replaced missing alleles based on mode versus mean showed that the structure and ranking of the top 9 models were identical, so we report only mode-based analyses.For the full dataset (i.e., all comparisons, without subsetting by distance), an additive model including the inverse of mean annual precipitation (1970-2000, a resistance model where aridity was linearly related to increased resistance) and the difference in average monthly maximum temperature Our non-invasive population, individual and individual-based landscape genetic analyses of aardvarks in southern and eastern Africa yielded insights on causes of genetic differentiation at scales ranging from local to continental.Across South Africa and Eswatini, arid areas increased resistance to gene flow.Differences in maximum temperatures at sampling sites also predicted individual genetic differentiation, suggesting the potential for isolation by adaptation, although that effect was apparent only after controlling for the effects of aridity and distance.While we did not evaluate climate change scenarios, these findings suggest that climate change will increase habitat fragmentation and limit gene flow for aardvarks, particularly where precipitation is expected to decrease and temperatures increase, on top of the adverse physiological and demographic consequences of drought observed by Rey et al. (2017) and Weyer et al. (2020).

F
Comparison of individualbased univariate landscape models of aardvark genetic variation in South Africa and Eswatini, based on AIC weight, for (A) all individual comparisons at geographic distances greater than lower thresholds ranging from 0 (all samples) to 1200 km and (B) all individual comparisons at geographic distances less than upper thresholds ranging from 20 to 1600 km (all samples).the mechanism further was difficult.Other factors identified in our univariate explorations of scale offered further support for the influence of climate on gene flow, such as the role of moderate temperatures at fine scales and moderate elevations at distances >200 km in facilitating gene flow.As noted by Van Buskirk and van Rensburg (2020), considering resistance versus difference measures appears informative as a signal of potential adaptive differences; difference measures would also yield useful alternative hypotheses when considering the consequences of recent fragmentation.Likewise, our use of distance thresholds to explicitly evaluate how the individual-based univariate landscape genetic models were influenced by scale provides Genetic differentiation among sampling areas separated by up to 430 km in central South Africa was relatively weak (Figures 2 and S2, individual-based; Tables

F
Genetic differentiation of individual aardvarks by distance, including (A) Mantel correlogram of genetic differentiation of individual aardvarks in South Africa and Eswatini, based on PCA distances estimated from 19 microsatellite loci, where significance (p < .05)and sign of correlation coefficients indicate more (positive) or less (negative) genetic similarity at those distance classes than expected; (B) average values of Lynch and Ritland (1999) relatedness as estimated for the combined central South African populations (dark triangles) and Kruger National Park (open circles), with values significantly (p < .05)different from 0 indicated with an asterisk (*).Sample sizes for each distance class were 39-42 (central South Africa) and 65-68 (Kruger); distances represent the maximum value of each distance class.
As aridity is predicted to increase in southernmost Africa under most climate change scenarios (The World Bank Group, 2022), we recommend more efforts to survey for aardvarks, identify areas where they have been extirpated or are absent, and establish baselines for distribution and monitoring, particularly in arid regions.Furthermore, the observed potential for adaptive differences associated with temperature underscore the critical need to resolve phylogeographic and subspecies questions across sub-Saharan Africa.Conservation assessments of the aardvark have, hitherto, been hampered by the lack of any modern evaluation of intra-species differentiation and adaptation to different habitats.Failing to recognize taxonomic divisions within the currently recognized species increases the risk of erroneously concluding conservation needs are low.

Locus name Forward primer Reverse primer Dye label Range (bp) Multi-plex Repeat motif NR a
a Number of repetitions in the reference genome.

Sampling area Country Samples collected Samples attempted to be genotyped Samples genotyped at ≥8 loci
Note: Count of samples genotyped successfully includes some individuals sampled multiple times, as ascertained by genotype matching.
. Difference predictors, including ΔT max in the top model isolation by adaptWeyer et al., 2020)th climate.Hot, dry and cold conditions cause physiological and behavioural changes in aardvarks(Rey et al., 2017;Weyer et al., 2020), suggesting the potential for selection around climate-related factors is high.As temperature and elevation difference predictors were correlated (TableS1), clarifying