Topographic barriers drive the pronounced genetic subdivision of a range‐limited fossorial rodent

Due to their limited dispersal ability, fossorial species with predominantly belowground activity usually show increased levels of population subdivision across relatively small spatial scales. This may be exacerbated in harsh mountain ecosystems, where landscape geomorphology limits species' dispersal ability and leads to small effective population sizes, making species relatively vulnerable to environmental change. To better understand the environmental drivers of species' population subdivision in remote mountain ecosystems, particularly in understudied high‐elevation systems in Africa, we studied the giant root‐rat (Tachyoryctes macrocephalus), a fossorial rodent confined to the afro‐alpine ecosystem of the Bale Mountains in Ethiopia. Using mitochondrial and low‐coverage nuclear genomes, we investigated 77 giant root‐rat individuals sampled from nine localities across its entire ~1000 km2 range. Our data revealed a distinct division into a northern and southern group, with no signs of gene flow, and higher nuclear genetic diversity in the south. Landscape genetic analyses of the mitochondrial and nuclear genomes indicated that population subdivision was driven by slope and elevation differences of up to 500 m across escarpments separating the north and south, potentially reinforced by glaciation of the south during the Late Pleistocene (~42,000–16,000 years ago). Despite this landscape‐scale subdivision between the north and south, weak geographic structuring of sampling localities within regions indicated gene flow across distances of at least 16 km at the local scale, suggesting high, aboveground mobility for relatively long distances. Our study highlights that despite the potential for local‐scale gene flow in fossorial species, topographic barriers can result in pronounced genetic subdivision. These factors can reduce genetic variability, which should be considered when developing conservation strategies.


| INTRODUC TI ON
The genetic subdivision and diversity of a species across space are determined by the combined effects of the environment and the ability of a species to disperse (Berthier et al., 2005;Manel et al., 2012;Quaglietta et al., 2013;Ruiz-Gonzalez et al., 2015).Dispersal ability can be limited by topographic barriers such as mountains or steep slopes, and by species-specific abiotic or biotic requirements, such as temperature or food availability, which may prevent the continuous distribution of individuals and reduce gene flow (Boulangeat et al., 2012;Cunningham et al., 2016;Cushman & Lewis, 2010;Sexton et al., 2009).As a consequence, natural selection, genetic drift, and inbreeding in smaller isolated populations may lead to heterogeneous patterns of genetic variability and population subdivision (Wright, 1969).Species with low dispersal ability such as those with low mobility and a burrowing lifestyle are especially prone to these processes.
Fossorial rodents build elaborate underground burrow systems.In many species, activities including searching for mates, reproduction, and foraging occur below-ground (Nevo, 1999).
Therefore, these rodents are often restricted to specific soil types and available food resources (Begall et al., 2007;Nevo, 1999;Reichman, 1975).Apart from these notable constraints in habitat and resource availability, the low mobility of fossorial species leads to small home ranges and limited dispersal (Harestad & Bunnel, 1979;Tucker et al., 2014).As a result, fossorial rodents often have a localized and patchy distribution.Moreover, in the case of solitary species, mature individuals meet mainly during the mating season, which further limits conspecific encounters.
These genetic and ecological patterns may be exacerbated in harsh environments, such as in mountain ecosystems, where the geomorphology of the landscape and the availability of suitable habitats limit dispersal opportunities and lead to restricted species distribution ranges and small effective population sizes (Badgley et al., 2017;Brown, 2001;Gaston, 2003;Rahbek, Borregaard, Colwell, et al., 2019).As a result, mountain regions have been recognized as hotspots for genetic differentiation and speciation, contributing disproportionately to terrestrial biodiversity, at least in the tropics (Rahbek, Borregaard, Antonelli, et al., 2019;Rahbek, Borregaard, Colwell, et al., 2019;Sandel et al., 2011).
However, species with limited distribution ranges and small population size, such as those found in mountain ecosystems, are at particular risk of extinction (Davies et al., 2009;Gaston, 2003).
Small populations tend to exhibit accumulations of deleterious mutations, low intraspecific diversity, and low adaptive potential, making them susceptible to environmental change and habitat shifts (Hoffmann et al., 2017;Lande, 1988;Willi et al., 2006).
Additionally, in mountain species, upslope range shifts in response to environmental change, and habitat shifts are limited by mountaintops (Parmesan, 2006;Wilson & Gutiérres, 2016).Mountain ecosystems face increasing threats from land use and climate change-induced habitat shifts.Therefore, it is imperative to understand the impact of species-environment interactions on genetic diversity in these ecosystems, especially for species with limited distribution and dispersal, and their implied threat of extinction.However, a thorough understanding is still lacking and this knowledge gap is particularly pronounced in African mountain ecosystems, where there is a scarcity of data from high-elevation research areas.Studying species in mountain ecosystems poses a significant challenge, especially for fossorial rodents, due to the difficulty of accessing remote areas and the inherent challenges in assessing these species in their natural below-ground habitats.
Our study addresses this knowledge gap by elucidating how landscape features drive the genetic subdivision and diversity of the giant root-rat (Tachyoryctes macrocephalus), a fossorial rodent endemic to the afro-alpine and afro-montane ecosystems of the Bale Mountains in southeast Ethiopia (Figure 1).
The species has a limited distribution range of ~1000 km 2 across the Bale Mountains massif and is found at altitudes ranging between 3000 and 4150 m above sea level (a.s.l.) (Sillero-Zubiri et al., 1995;Yalden & Largen, 1992).Giant root-rats have specific habitat requirements, occurring in grasslands in areas with good soil depth, especially along wetland shores and flooded valleys (Sillero-Zubiri et al., 1995;Šklíba et al., 2017).Grasslands in river valleys that spread through shrubs and forest zones into lower elevations allow the species to expand down to about 3000 m a.s.l.(Yalden, 1985).Their relatively small home ranges (about 100 m 2 ) can shift throughout the year depending on food availability (Šklíba et al., 2020).Giant root-rats are significant ecosystem engineers and create large underground burrow systems, in which they live solitarily (Yalden, 1985).Through their combined effect of soil perturbation and herbivory, they alter nutrient availability, soil texture and moisture, and maintain their own habitat by keeping the vegetation low (Asefa et al., 2022;Miehe & Miehe, 1994;Šklíba et al., 2017;Yalden, 1985).By using below-ground burrows, the species circumvents the harsh environmental conditions of the mountain ecosystem, which include strong winds and temperatures below 0°C, and limits the risk of being preyed upon by its main predator -the Ethiopian wolf (Canis simensis) (Sillero-Zubiri & Gottelli, 1995;Šumbera et al., 2020;Vlasatá et al., 2017; afro-alpine, conservation genetics, fossorial rodents, landscape genetics, population genetics Yalden, 1985).This combination of the species' predominant below-ground activity, limited dispersal and distribution range within a changing African mountain ecosystem, positions the giant root-rat as an ideal model organism to investigate the connection between genetic patterns and landscape features.This exploration aims to better understand how environmental change impacts population subdivision in remote mountain ecosystems for species with limited dispersal, so as to preserve mountain biodiversity and ensure ecosystem functioning.
In the present study, we analysed spatial genetic subdivision and diversity in the giant root-rat across its distribution range.To achieve this, we analysed both mitochondrial genomes (mitogenomes) and nuclear genomes and further investigated the relationship between genetic differentiation and landscape features.
We generated complete mitogenomes and low-coverage nuclear genomes from 77 individuals collected across nine sampling localities in the Bale Mountains (Figure 1c).We applied two different landscape genetic approaches to evaluate how gene flow of the species is impacted by geographic distance, by vegetation and soil moisture, and by slope and elevation (used as proxies for topographic barriers).Due to the predominantly below-ground activity and patchy distribution of the giant root-rat, we hypothesize strong genetic subdivision across small spatial scales.Owing to the pronounced heterogeneity of the environment in the Bale Mountains, we hypothesize that genetic structuring is driven by habitat availability and by topographic structures across the species' range.

| Study area
The Bale Mountains in southeast Ethiopia (6°29′ N-7°10′ N and 39°28′ E-39°57′) represent Africa's largest afro-alpine ecosystem, comprising ~8% of the continent's area above 3000 m a.s.l.(Groos et al., 2021, Figure 1a-c).In order to protect the unique afro-montane and afro-alpine ecosystems of the Bale Mountains, the area above ~3200 m a.s.l. was designated a national park in 1970.The Bale Mountains are characterized by a dry season from November to February, and a rainy season with two parts consisting of short rains from March to June and long rains from July to October (Miehe & Miehe, 1994).The vegetation of the Bale Mountains shows an elevational zonation from moist montane forest (~1500-3500 m a.s.l.) over ericaceous shrubland and dwarf forest (~3500-4000 m a.s.l.) to afro-alpine vegetation with open grassland and Erica outposts (above 4000 m a.s.l.).Dwarf-shrub vegetation, such as Helichrysum associated with Lobelia, is the main plant formation in the afro-alpine vegetation but does not cover the whole area, leaving open spaces for herbaceous plants like Senecio, Alchemilla, or Salvia (Miehe & Miehe, 1994;Tallents & Macdonald, 2011).Late Pleistocene, between 42,000 to 16,000 years ago (kya, Groos et al., 2021).The plateau is bounded by several outlet valleys in the north and the east with slopes covered by dense, shrubby Erica vegetation and by congealed lava flows at its northwestern margins (Figure 1c).These topographic structures distinguish the plateau from the northern region of the national park, which is ~300-500 m lower in elevation and comprises broad valleys and plains with afro-alpine vegetation (Miehe & Miehe, 1994).In comparison to the plateau, the north has higher moisture availability and milder temperatures.

| Sampling
We collected tissue samples from 77 live giant root-rat individuals at nine localities across the Bale Mountains National Park, covering the entire distribution range of the species (Figure 1c, Table S1).The sampling localities were distributed across the two topographically distinct regions in the national park, in the north (localities N1-N3) and in the south (localities S1-S6).The southern localities are scattered across the central and southern parts of the Sanetti Plateau.Localities sampled in the north of the plateau lie at a lower elevation (~3500 m a.s.l.) than localities sampled in the south (~3800-4000 m a.s.l.).Sampling localities between regions were separated by 15.3-28.3km, and localities within regions were separated by 2.6-16.0km.We captured 7-9 giant root-rat individuals per locality (except for locality S4 with n = 1).The samples were collected in January and February in two consecutive years (2020: n = 60, 2021: n = 17) under the permit of the Ethiopian Wildlife Conservation Authority.Individuals were caught with snare traps that were monitored by the capture team at all times to guarantee no harm to the animals.An ~0.5 cm 2 piece of skin from the hind leg was cut with sterilized scissors and stored in 96% ethanol or DNAgard® for blood and tissue (Biomatrica, Inc.) for genomic analyses.After sterilizing the wound, the animals were immediately released back into their burrow systems.

| Laboratory analyses
We extracted DNA from the tissue samples using the Qiagen DNeasy® Blood and Tissue Kit following the manufacturer's protocol (Qiagen Ltd.).Sixty of the samples collected in the 2020 field season were processed in-house in the modern DNA labs at Globe Institute, University of Copenhagen.The DNA concentration of the extracts was measured using QubitTM dsDNA HS (Invitrogen).After quantification, we diluted the extracts to a concentration of 6 ng/μL in a total volume of 50 μL.
DNA was sheared to ~400 base pair (bp) fragment lengths using the Covaris M220 ultrasonicator.We built DNA fragments into an Illumina library following the protocol from Carøe et al. (2018) and doubleindexed them using AmpliTaq Gold Polymerase (ThermoFisher) during the indexing PCR step.Index PCR reactions were performed in 100 μL, using 1× PCR buffer, 2.5 mM of MgCl 2 , 0.2 mM of dNTPs, 0.2 μM of index primer mix, 0.1 U/μL of polymerase, and 20 μL of DNA library.PCR cycling conditions were 95°C for 10 min; 10-18 cycles of 95°C for 30 s, 50°C for 30 s, and 72°C for 1 min followed by 72°C for 5 min.
The number of cycles for the index PCR was determined from qPCR analysis, using the threshold cycle (or Ct value) as the metric to determine the number of cycles.Post-amplification, libraries were purified using SPRI beads as in Carøe et al. (2018).The purified indexed libraries were quantified on a QubitTM dsDNA HS (Invitrogen), and qualitychecked on either an Agilent 2100 Bioanalyzer or an Agilent Fragment Analyzer™.Libraries were pooled equimolarly and sequenced on an Illumina NovaSeq 6000 using paired-end (PE) 150 bp technology (Novogene Europe, http:// en.novog ene.com).
For the remaining 17 samples collected in the 2021 field season (localities S1 and S2), DNA was extracted and processed to libraries by Novogene to expedite sample processing, and sequenced on a NovaSeq 6000, using 150 bp PE technology.

| Data generation of mitochondrial and nuclear DNA
We trimmed adapters and removed reads shorter than 30 bp for each individual using skewer v0.2.2 (Jiang et al., 2014).We merged overlapping paired-end reads using FLASHv1.2.11 (Magoč & Salzberg, 2011) with default parameters.We mapped both merged and unmerged reads to the hoary bamboo rat (Rhizomys pruinosus) nuclear genome (Genbank accession: VZQC00000000.1;Guo et al., 2021) which is the nearest relative of the giant root-rat with an available genome, combined with the giant root-rat mitogenome (Genbank accession: MW751806; Reuber et al., 2021).The divergence time between the hoary bamboo rat and the giant root-rat is ~14 million years (median time estimate based on 4 studies from www. timet ree.org), and mapping to a divergent reference may lead to the recovery of more conserved regions in the genome.However, as all of our individuals are a single species and equally distant to the hoary bamboo rat, they should all be biased to the same extent, meaning our overall inferences should not be impacted.We used BWA v0.7.15 (Li & Durbin, 2009) utilizing the mem algorithm and default parameters.We parsed the alignment files, and removed duplicates and reads of mapping quality score <30 using SAMtools v1.6 (Li et al., 2009).We built consensus mitogenomes from each individual using a majority rules approach (-doFasta 2) in ANGSD v0.921 (Korneliussen et al., 2014) only considering bases with a base quality score greater than 30 (-minq 30), reads with a mapping quality score greater than 30 (-minmapq 30), and sites with at least 10× coverage (-minInddepth 10).The mitogenomes are available under GenBank accessions OQ207545-OQ207620.

Phylogenetic analysis
We used PartitionFinder v. 2.1.1 (Lanfear et al., 2017) to find the best model of evolution for the mitochondrial genome data.The Bayesian phylogeny based on the 48 mitogenome haplotypes identified in the network analysis was constructed using the GTR + I + G model of evolution and MrBayes v.3.2.7a (Ronquist & Huelsenbeck, 2003).The MCMC algorithm was run twice with four chains of 10 million generations, sampled every 1000 generations and with a 10% burn-in.The trees were combined following the majority-rule consensus approach, to assess the posterior probability of each clade.The resulting tree was midpoint rooted and visualized in FigTree v.1.4.4 (http:// tree.bio.ed.ac.uk/ softw are/ figtr ee/ ; Figure S1).

Fixation statistics and AMOVA
We calculated the pairwise differentiation between eight of the nine sampled localities (omitting S4 as n = 1) with the F ST -estimator of the software Arlequin v. 3.5.2.2 with 10,000 permutations (Excoffier & Lischer, 2010); p-values of F ST -estimates were adjusted using the Bonferroni correction (Rice, 1989), controlling for a false-positive discovery rate (R Core Team, 2021).Additionally, we conducted a hierarchical analysis of molecular variance (AMOVA), also in Arlequin, with localities grouped into their provenance in the regions north and south (Figure 1c).

| Nuclear DNA analysis
We investigated population subdivision of the nuclear data using principal component analysis (PCA) and admixture proportion analysis.We generated genotype likelihoods in ANGSD (Korneliussen et al., 2014) for all individuals using the following filters and parameters: call genotype likelihoods using the GATK algorithm (-GL 2), output a beagle file (-doGlf 2), only include reads with mapping and base qualities greater than 30 (-minmapQ 30 and -minQ 30), only include reads that map to one location uniquely (-uniqueonly 1), a minimum minor allele frequency of 0.05 (-minmaf 0.05), only call an SNP if the p-value is greater than 1e −6 (-SNP_pval 1e-6), infer major and minor alleles from genotype likelihoods (-doMajorMinor 1), only include sites if at least 40 individuals are covered (-MinInd 40), remove scaffolds shorter than 100 kb (-rf), and remove secondary alignments (-remove_bads 1).To compute the PCA, we constructed a covariance matrix from the genotype likelihoods using PCAngsd v0.98 (Meisner & Albrechtsen, 2018).Admixture proportions were calculated using the same genotype likelihoods with NGSadmix (Skotte et al., 2013).
We ran NGSadmix specifying K = 2 and K = 3.To evaluate the reliability of the NGSadmix results, we ran each K up to 100 times independently.If we retrieved consistent log-likelihoods from at least two independent runs, the corresponding K was considered reliable.
To estimate levels of differentiation among localities, we computed F ST -estimates from a consensus haploid call file created using ANGSD (-dohaplocall 2) and the same filtering parameters as the PCA and admixture proportions above.We calculated the F ST -estimates using an available python script (https:// github.com/ simon hmart in/ genom ics_ gener al/ blob/ master/ popge nWind ows.py) and specifying a window size of 1 Mb, and a minimum number of sites per window as 1000 bp.

Ancient gene flow between regions
The mitogenomes of two individuals (WM07 from locality S1 and GG01 from locality S6) grouped with individuals from the north in the haplotype network and phylogeny.We therefore used Dstatistics (also known as ABBA/BABA, Durand et al., 2011) to test whether the results were driven by ancient gene flow between regions north and south, or by incomplete lineage sorting.We tested several topologies [[H1, H2], H3], with branch H1 being one of the two putative introgressed individuals, and branches H2 and H3 being individuals from one region north or south, or one from each region.A negative D-score illustrates a closer relationship between H1 and H3 than H2 and H3, while a positive D-score indicates that branches H2 and H3 are more closely related than H1 and H2 (for an example see Table 1).This setup can also be used to uncover population subdivision, as the incorrect input topologies would lead to elevated D-scores due to more recent common ancestry, as opposed to gene flow (Westbury et al., 2018).We performed the D-statistic tests using a random base-call approach in ANGSD (-doAbbababa 1) (Korneliussen et al., 2014).We implemented the same filtering approach as for the above analyses but only included scaffolds >1 megabase (Mb) in size, a block size of 1 Mb (-blocksize 1000000), and the hoary bamboo rat (Rhizomys pruinosus, GenBank accession VZQC00000000.1,Guo et al., 2021) as the ancestral state/outgroup (-anc).To assess the significance of our results we used a block jackknife test with the script jackKnife.R which is available with the ANGSD tool suite.

| Diversity
Based on the mitochondrial genomes, we calculated nucleotide diversity (π) per region (north and south), and separately for eight of the nine localities using DnaSp v.6 (omitting S4 as n = 1) (Rozas et al., 2003).We tested differences in nucleotide diversity between the two geographic regions, and among localities using genetic_di-versity_diffs v. 1.0.3 (Alexander, 2017).
The python script used to calculate nuclear F ST above simultaneously computes nucleotide diversity per region and per locality.
To test for significant differences in levels of nucleotide diversity between regions and between localities, we used a Welch-test (unpaired t-test), accounting for unequal variance.

| Landscape genetic analysis
We applied two landscape genetic approaches to investigate the effects of landscape features on the observed genetic differentiation between localities based on F ST -estimates of the mitochondrial and nuclear genome, using partial Mantel tests and resistance surface optimization.The limited number of sampling localities prevented us from analysing the north and south regions separately.
We selected four environmental variables (vegetation, soil moisture, slope, elevation), which were based on satellite-based remote sensing data, as predictors for genetic differentiation of the giant root-rat.For vegetation and soil moisture, we used annual mean raster layers of the Normalized Differentiation Vegetation Index (vegetation index; Bannari et al., 1995) and the Land-Surface Water Index (soil moisture index; Chandrasekar et al., 2010;Gao, 1996).These raster layers were calculated based on Sentinel-2 satellite imagery, which was captured twice per month throughout 2021, downloaded from ESA's Copernicus Hub (https:// scihub.coper nicus.eu/ ), preprocessed using the FORCE workflow (David Frantz, 2019) and provided by Bhandari et al. (2023).The near-infrared (band 08A) and shortwave infrared (band 11) satellite bands were used to calculate the soil moisture index (Chandrasekar et al., 2010;Gao, 1996).The In that way, we were able to account for high correlation among ma- for the non-independence of the predictor variables and to account for spatial autocorrelation (Clarke et al., 2002;Peterman et al., 2014;Shirk et al., 2018).This iterative process works towards identifying the best-fit landscape resistance surface.
In our optimization framework, we used a single surface optimization approach, where the resistance surfaces of the generated raster layers of each selected environmental variable (i.e.vegetation index, soil moisture index, slope, and elevation, see paragraph "Landscape genetic analysis") were optimized individually, using eight transformation functions (Monomolecular and Ricker functions) and the default parameters (Peterman, 2018).In this step, the pairwise resistance distance between localities was estimated by assuming that individuals can use several paths to disperse.Resistance distances were generated with the costDist function implemented in the ResistanceGA package and movements between localities were allowed in eight directions during resistance distance calculation.Because Euclidean distance is incorporated in the resistance distances, it was not included as an additional variable.In the optimization process, the mixed-effect models were calculated, fitting the pairwise genetic differentiation as a response against the resistance distances as single fixed effects, using the AIC for model evaluation and including sampling localities as a random effect to account for spatial autocorrelation (Table 2).
Following the single surface optimization approach, we conducted multisurface resistance optimizations of non-correlated resistance surfaces, following the workflow of Peterman (2018).Therewith, we TA B L E 2 Model selection results for linear mixed-effect models, testing the effect of resistance distances (dispersal costs based on environmental condition) on levels of (A) mitochondrial and (B) nuclear genetic differentiation (F ST ) in the giant root-rat.

| RE SULTS
We generated complete mitogenomes from all 77 giant root-rat individuals, with a depth of coverage ranging between 25× and 383×, and a length of 16,646 bp.For the nuclear genomes, we obtained coverages ranging from 0.12× to 0.77× (Table S2).SNPs calling resulted in 8,557,525 sites.

| Haplotype network and phylogeny
Our 77 sampled giant root-rat individuals comprised 48 haplotypes, with no haplotypes shared among localities (Figure 2a).The haplotype network and phylogeny of the mitogenomes revealed two geographically separated and well-supported groups; one group comprising all samples collected in the north, with the inclusion of individuals WM07 (locality S1) and GG01 (locality S6) from the south, and one group comprising all other samples from the south (Figure 2a; Figure S1).In the network, we identified 574 segregating sites among individuals, and the northern and southern haplogroups were separated by 334 segregating sites (Figure 2a).The north contained 32 individuals (21 haplotypes) and the south contained 45 individuals (27 haplotypes).Within the north, we identified two distinct, well-supported genetic clades, N A and N B (Figure S1).Clade N A comprised six individuals (five haplotypes) from localities N1-N3.
Individual GG01 from locality S6 was basal to the clade.Clade N B comprised the remaining individuals from the north with individual WM07 from the south at basal position.We did not identify any spatial genetic structuring in the south.
The AMOVA yielded a high level of between-region variation when the eight localities (omitting S4 with n = 1) were grouped into north and south (89.50%, p < .01).Within each region, the variation was higher within localities (9.27%) than among localities (1.24%).

| Principal component analysis and admixture proportions
We identified two main groups in the nuclear data, in agreement with the mitochondrial findings (Figure 2).In the PCA, individuals from the north separated from the individuals from the south, with almost 10% of the variation explained on the first principal component.The southern group showed a slight separation on the second component, with the more central localities S1 and S2 segregating from the localities further southeast, with 1.4% variation explained (Figure 2c).The division of the data into north and south was also evident in the admixture analysis of K = 2 (Figure 2d).There was no sign of recent admixture between the two southern individuals WM07 and GG01 and the north, contradicting the mitochondrial lineage pattern.The admixture analysis did not converge with K = 3, suggesting K = 3 did not reliably fit the data.

| Ancient gene flow between regions
To investigate the origin of the mitochondrial lineages, present in individuals WM07 and GG01, which were more closely related to the northern haplogroup than the south (Figure 2a), we tested for ancient gene flow using the nuclear data.Using the topology [[WM07,south], north], we found most comparisons to have a D-score around 0 and a Z-score < |3|, indicating that WM07 and all individuals from the south were equally related to the northern individuals (Figure 3, Figure S2).We found positive D-scores (Z-score >3) when using the topology [[WM07,north], north], demonstrating a closer relationship between individuals from the north with each other than with WM07, which agrees with their more recent common ancestry and the basal position of WM07 in the phylogenetic tree (Figure S1).We found qualitatively the same results when we investigated the relationship of GG01 (sampled in locality S6) with individuals from the north and south (Figure 3, Figure S2).Hence our analysis did not support the mitochondrial lineages in WM07 and GG01 as the result of ancient gene flow between north and south.

| Genetic differentiation among localities
Investigating pairwise genetic differentiation between localities, we found that F ST -estimates were higher between regions north and south, than within regions, at both the mitogenome and nuclear level (Figure 4).For the mitogenomes, pairwise differences between localities of different regions ranged from 0.82 to 0.97 and were much higher than within regions, where values ranged from 0.09 to 0.29.
For the nuclear data, pairwise F ST -estimates between localities from different regions ranged from 0.04 to 0.06 and were higher than values between localities within regions, which ranged from 0.01 to 0.02 (Figure 4).

| Diversity
We estimated levels of diversity for each region, and for each locality (omitting S4 as n = 1).For the mitogenomes, diversity in the south (π = 0.003 ± 0.0002) was significantly higher than in the north (π = 0.001 ± 0.001, p < .05; Figure 2b), which reflected the presence of the divergent mitochondrial lineages in individuals WM07 and GG01 from localities S1 and S6 (Figure 2a).Thus, this was also apparent in localities S1 and S6 having the highest diversity levels, with S6 showing significantly differentiated levels of diversity (p < .05)and S1 showing marginal significant differentiation (p > .05< .1,Table S3) to the remaining localities in the south.
Based on the nuclear data, we also observed significantly higher diversity in the south (π = 0.257 ± 0.0364) than in the north (π = 0.249 ± 0.0357, p < .05; Figure 2b).Among localities in the south, the two central localities S1 and S2 had lower diversity levels than the localities sampled further to the southeast S5, S6).All localities had significantly differentiated levels of diversity (p < .05,Table S3).

| Landscape genetics
We used a reciprocal causal modelling approach with partial Mantel tests to relate genetic differentiation to the distance matrices of the environmental variables vegetation index, soil moisture index, slope, and elevation.For mitochondrial and nuclear genomes, we found that elevation was the strongest model.Elevation was significantly related to genetic differentiation in all partial Mantel tests, regardless of the second environmental distance matrix controlling elevation (Tables S4A, S5A).Furthermore, elevation showed the strongest relative support; all relative correlation values were positive, after the effect of the other environmental models was removed (Tables S4B, S5B).Geographic distance was the second strongest model and showed a significant relation with genetic differentiation when controlled by vegetation index, soil moisture index, or slope.
For mitochondrial genomes, geographic distance was non-significant when controlled by elevation, while elevation was significant with a negative correlation value at the nuclear genome level.Slope was significantly correlated with mitochondrial genetic differentiation when controlled by geographic distance, vegetation index, and soil moisture index, but showed a negative correlation value when controlled by distance, and a non-significant correlation value when controlled by elevation (Table S4).The nuclear genomes revealed similar results (Tables S5).
For mitochondrial and nuclear genomes, we could not find significant effects of vegetation and soil moisture indices on genetic differentiation; all variables had non-significant Mantel correlation values, showing the least support in the reciprocal causal modelling matrix (Tables S4, S5).
The raster layer optimization approach revealed that slope and elevation were the best-fitting models for explaining genetic differentiation at both mitochondrial and nuclear levels.For the mitogenomes, slope was selected as top model in 69.4%, elevation in 27.3%, and geographic distance in 3.3% across 1000 bootstrap iterations (Table 2a).For nuclear genomes, the top model was slope in 60.7% of the iterations and elevation in 39.3% (Table 2b).The average weight and the AICc from the bootstrap analyses supported geographic distance as a driver for genetic differentiation, while the parameters rank, AIC, maximum likelihood, and R 2 suggested slope and elevation as the best models explaining mitochondrial and nuclear genome genetic differentiation (Table 2).The response curve of the optimized resistance surface showed that the resistance costs increased with increasing slope (Figures S3d, S5d

| Genetic subdivision between regions
The significant north and south geographic subdivision in the mi- The presence of genetically distinct subpopulations may be attributed to long-lasting extrinsic barriers, which prevent genetic exchange between them (Avise, 2000;Bryja et al., 2010).Slope and elevation were identified as the primary drivers of genetic differentiation between north and south in our landscape analysis (Table 2, Tables S4, S5), similar to what has been observed in other fossorial rodents such as the Brazilian tuco tuco of the dunes (Ctenomys flamarioni), and the common water vole (Arvicola terrestris) (Berthier et al., 2005;Fernández-Stolz et al., 2007).The south comprises the Sanetti Plateau, which is ~300-500 m higher in altitude than the northern region (Figure 1c).The plateau margins northwards are characterized by broad valleys with steep slopes covered by dense Erica thickets.In addition to the slopes, which themselves act as a barrier, the Erica thickets may further limit the dispersal of the species (Miehe & Miehe, 1994;Yalden, 1985).Giant root-rats are adapted to open grasslands.They avoid dense shrubs such as Erica (Yalden, 1985), likely due to the difficulties of burrowing in woody ground and the absence of food-plants.Additionally, the plateau is bounded by congealed lava flows of unknown age to the northwest.
These barriers to dispersal and the fossorial lifestyle of the giant root-rat, which limits the species' ability to traverse pronounced topographic structures, presumably caused the strong genetic subdivision of the species.Our landscape genetic findings align with recent satellite-based predictions of giant-root rat distribution, which indicated the species' distribution can be explained by texture metrics of the landscape, i.e. factors that characterize topographic variation across the landscape (Wraase et al., 2022).
In addition to slope and elevation, the pronounced subdivision observed in giant root-rats may also have been reinforced by glacial extents during the Late Pleistocene.The Bale Mountains are currently ice-free, but the Sanetti Plateau in the south was glaciated between ~42 to 16 kya (Figure 1c;  S3).This would be in agreement with the often-proposed hypothesis that populations of mammals exhibit reduced genetic diversity on recently deglaciated land (e.g.Hewitt, 1996Hewitt, , 2004)).The glacial extent in the north and northwestern valleys of the plateau margins persisted until ~16 kya, while the ice shield on the plateau around Tullu Dimtu was smaller in extent already ~20 kya (Groos et al., 2021).
Although vegetation and soil moisture were previously identified as essential factors influencing the local abundance of giant root-rats (Asefa et al., 2022;Šklíba et al., 2017), our study did not indicate any effect on genetic differentiation (Table 2).This suggests these factors play less of a role in hindering gene flow at the landscape scale between regions.However, the spatially coarse vegetation and soil moisture indices used in our analysis may not fully capture the highly specific food and soil requirements of the giant root-rat.Their primary food resource is Alchemilla (Yaba et al., 2011).The vegetation index, which is based on remotely-derived satellite data, may not distinguish its spectral signal from other preferred plants (Wraase et al., 2022).Additionally, the giant root-rat requires soil layers of approximately 50 cm in depth to engineer burrow systems and for thermoregulation (Sillero-Zubiri et al., 1995;Šumbera et al., 2020).While

| Gene flow within regions
We observed a lack of structuring among localities within the regions north and south at the local scale.Levels of differentiation were low, with nuclear F ST -estimates of 0.01-0.02within regions, which is considered weak differentiation for nuclear data (Figure 4; Weir & Cockerham, 1984;Wright, 1978).This indicates a high level of dispersal and gene flow across distances of at least 16 km, which was the maximum distance between two sampling localities within regions.The ability of giant root-rats to disperse across such relatively large distances was in contrast to our expectations; giant root-rats are fossorial, solitary, and territorial (Yalden, 1985).In combination with the heterogeneity in soil structure and food availability across its range, we had expected this would lead to stronger genetic structuring at small spatial scales, similar to what has been reported for other fossorial rodents (Mapelli et al., 2012;Schweizer et al., 2007).
Although direct observations for this are still lacking, the limited substructuring within regions at the local scale, and the large dispersal distances suggest giant root-rats can disperse aboveground and across relatively large distances.In fact, giant root-rats show morphological adaptations to surface activity, in that their eyes are situated dorsally on the head (Figure 1d), which allows them to detect predators in open habitats (Yalden, 1985).In support of our findings, radio tracking has evidenced the dispersal of a giant root-rat individual over a distance of up to 270 m within a span of 2 days; the tracked individual traversed across damp soil, suggesting it did not disperse underground (Šklíba et al., 2020).Aboveground dispersal has also been documented in other fossorial, solitary rodent species, such as blind mole-rats (Spalax microphthalmus; Zagorodniuk et al., 2018) and Tibetian plateau zokors (Eospalax fontanierii; Chu et al., 2021).Even in strictly subterranean African mole-rats, long-distance dispersal is not precluded (Fukomys damarensis, Bathyergidae; Finn et al., 2022).
For giant root-rats, aboveground dispersal attempts could be triggered by decreasing food supply, the absence of sexual partners, or the presence of competitors (Šklíba et al., 2020;Zagorodniuk et al., 2018).Also, the behaviour may circumvent the patchy availability of suitable habitats and small home-ranges, maintaining gene flow and limiting genetic structuring across small spatial scales.
Dispersal events in the giant root-rat may be male-dominated, as it has been observed in tuco tucos (Ctenomys talarum and C. australis; Cutrera et al., 2005;Mora et al., 2010), Chinese zokor (Eospalax fontanierii; Zhang, 2007), giant mole-rats (F.mechowii; Kawalika & Burda, 2007), and arvicoline rodents (Le Galliard et al., 2012).While sex-specific dispersal has not been studied in the giant root-rat yet, the observation of males being more frequently involved in dispersal attempts compared to females (Šklíba et al., 2020) combined with microsatellite analyses that indicate males disperse for longer distances suggests that this type of dispersal may be prevalent in the species.
Our nuclear data did suggest slight subdivision in the south, with localities in the central part of the plateau (S1, S2) being more differentiated from localities in the southeast (S3-S6; Figure 1c).This was evidenced by increased F ST in their pairwise comparisons and their segregation on the second principal component on the PCA (Figure 2c,d).Although overall differentiation among localities in the south was low, this pattern may reflect topographic features; the mountain Tullu Dimtu, the highest peak in the Bale Mountains National Park with 4377 m a.s.l., is located close to localities S3 and S4, and may limit gene flow (Figure 1c).

| Conservation implications
Through landscape genetic analysis, we identified the drivers of population subdivision between north and south to be topographic barriers in the form of slope and elevation.While the species is capable of dispersing locally, our findings suggest that giant root-rats in the north and south must be considered separately when developing conservation strategies, as we did not find evidence of landscape-scale gene flow between the two regions.Giant root-rats impact their surrounding environment as ecosystem engineers and are primary prey for the endangered Ethiopian wolf, which underscores the importance of their persistence (Sillero-Zubiri & Gottelli, 1995;Šklíba et al., 2017).Already, the giant root-rat is believed to have a small census size due to its limited distribution range (although no census estimate is available) and is listed as endangered by the IUCN (Lavrenchenko & Kennerley, 2016).
The potential for reduction of the species' distribution range due to increasing human activities in the form of expanding livestock grazing and human settlements in the Bale Mountains could harm the species' persistence with negative consequences for the overall ecological balance in the region (Gashaw, 2015;Mekonen, 2020;Stephens et al., 2001).Our study yields some key insights for planning future conservation strategies for the species and highlights the value of genomic data in expanding our understanding of the population dy- We are grateful to the Ethiopian Wildlife Conservation Authority Bale Mountains National Park is the afroalpine Sanetti Plateau, which spans elevations from approximately 3800 m a.s.l. to 4377 m a.s.l. at the peak of the mountain Tullu Dimtu (Figure 1c).Large parts of the plateau were glaciated during the F I G U R E 1 Sampling localities of giant root-rats within Bale Mountains National Park, Ethiopia.(a) Map of Africa showing Ethiopia in dark grey; (b) map of Ethiopia indicating the location of Bale Mountains National Park; (c) map of the nine sampling localities from two distinct geographic regions.The north (~3500 m above sea level [m a.s.l.]), and south (~3800-4000 m a.s.l., Sanetti Plateau) are separated by steep slopes covered with dense Erica thickets and congealed lava flows.Sample size of each locality is indicated in brackets.Tullu Dimtu is the highest peak in the Bale Mountains National Park at 4377 m a.s.l., and is indicated with a black triangle.Region with light blue shading indicates the glacial extent within Bale Mountains National Park ~35 ± 7.1 thousand years ago (kya) (Groos et al., 2021; Ossendorf et al., 2019); (d) burrowing giant root-rat, photography by V. Reuber.
vegetation index was calculated directly within the FORCE workflow and provided in theBhandari et al. (2023) dataset(Bannari et al., 1995).The vegetation and soil moisture index were used as proxies for food and soil availability for giant root-rats(Sillero-  Zubiri et al., 1995;Šklíba et al., 2017;Yaba et al., 2011;Yalden & Largen, 1992).To determine whether topographic structures act as barriers for burrowing giant root-rats, we included the variables slope and elevation in our analysis.Raster layers for slope and elevation were obtained from a Shuttle-Radar-Topography-Mission digital elevation model from the USGS Earth explorer (www.earth explo rer.usgs.gov).The generated raster layers of all four environmental variables had a 30 × 30 m resolution and were cropped to the extent 6°40′ N-7°02′ N, to 39°36′ E-39°57′ E. Our analyses were conducted in R environment version 4.2.1 (R Core Team, 2021).2.7.1 | Partial Mantel testsUsing partial Mantel tests, we analysed if the genetic differentiation of mitochondrial and nuclear genomes between eight of the nine sampling localities (omitting S4 as n = 1) was correlated with geographic distance, vegetation index, soil moisture index, slope, and elevation.Therefore, we constructed distance matrices, using the generated raster layers of the four environmental variables (see paragraph "Landscape genetic analysis").The genetic distance matrix was generated by linearizing the pairwise genetic differentiation estimates between localities, that is, the F ST -estimates ([F ST /(1-F ST )];Rousset, 1997).The geographic distance matrix was calculated in Euclidean distances and log-transformed to linearize the relationship with genetic distance.For each environmental variable, we extracted their values from the computed raster layers at the coordinates of the sampling localities and therewith generated the environmental distance matrices.We then applied a pairwise reciprocal causal modelling approach.Reciprocal causal modelling compares partial Mantel tests of a focal environmental model, removing the influence of a competing, alternative model(Cushman et al., 2006;Cushman & Landguth, 2010).In this approach, the correlation between one environmental distance matrix and genetic distance is controlled by a second matrix (e.g.focal model: genetic distance~geographic distance|elevation distance), and in a next step, both environmental distance matrices are interchanged (e.g.alternative model: genetic distance~elevation distance|geographic distance).
trices(Cushman et al., 2006;Cushman & Landguth, 2010).To assess which of the two models explains genetic distance better, the relative support of the focal and alternative model was calculated by estimating the difference between the correlation values of the two models.If the difference in correlation factors was positive, we assumed that the focal hypothesis was correct.The partial Mantel tests were performed with 9999 permutations in the vegan R package v.2.6-4(Oksanen et al., 2022).The R scripts and data files for the partial Mantel tests are available at: https:// zenodo.org/ recor ds/ 10401840.2.7.2 | Raster layer optimization framework to generate resistance surfacesWe used a raster layer optimization framework developed byPeterman (2018) andPeterman et al. (2014), to further identify landscape features that explain mitochondrial and nuclear genetic differentiation, using the R package ResistanceGA(Peterman, 2018).In this framework, the generated raster layers of the environmental variables (i.e.vegetation index, soil moisture index, slope and elevation) were transformed into resistance surfaces and optimized with the ResistanceGA package utilizing a genetic algorithm from the GA R package(Scrucca, 2013).A resistance surface is a spatial layer that assigns values to each grid in the raster layer of the selected environmental variable.Those values are used to estimate the cost of dispersal and mirror to what extent the selected variable hinders or facilitates the connectivity of a species between two localities (pairwise resistance distances).Thereby, there are no a priori assumptions about the relationship between the environmental variable and the species' dispersal characteristics.The genetic algorithm in the optimization framework is used to maximize the relationship between the resistance distances of each raster layer, and the pairwise genetic differentiation (F ST ) between localities.The process of generating resistance surfaces is repeated, and in every iteration, the resistance distances are fitted against the genetic distances in a mixed effect model, until the objective function, the AIC (Akaike's information criterion;Akaike, 1974) of the mixed effect model does not improve further.The mixed-effect models are conducted using a maximum likelihood population effect parameterization to account tested if several surfaces better represent the environment causing genetic differentiation of giant root-rats.The resistance surfaces of the vegetation index and soil moisture index were both correlated to elevation and with each other.Consequently, we performed three distinct multisurface optimization runs.In each run, we separately optimized the pairs of resistance surfaces slope and vegetation index, slope and soil moisture index, and elevation and slope.During the optimization process, each surface in the multisurface analysis was modified using the methods described earlier.The modified surface pairs were then combined to create a single, composite resistance surface.This composite surface was subsequently utilized to calculate pairwise effective distances.We combined all cost distance matrices of single and multisurface optimization outputs in a list, on which we conducted the bootstrap model selection with 75% of the samples and 10,000 iterations.We used the AIC as the objective model function during optimization and the AICc as the rank variable.The bootstrap model selection involves re-fitting the mixed-effect models and computing fit statistics for each model.This process reveals the average AIC and the percentage by which each resistance surface is chosen as the topranked model across all bootstrap(Peterman, 2018).The R scripts and data files for the raster layer optimization framework are available at: https:// zenodo.org/ recor ds/ 10401840.

F
Patterns of genetic diversity and subdivision of giant root-rats across their range.(a) Network of the 48 mitochondrial haplotypes present among the 77 sampled individuals.Each circle represents a haplotype, and the relative size of each circle represents haplotype frequency.Numbers on branches show the number of segregating sites between haplotypes for >10.Black dots indicate intermediate haplotypes not present in the data.Distances between haplotypes are not to scale.(b) Diversity levels within the sampled localities and subpopulations based on mitogenomes (circles) and nuclear data (triangles).Sample sizes are shown in brackets.(c) PCA based on 77 low-coverage nuclear genomes.The percentage of the diversity of each principal component is shown in brackets.(d) Ancestry proportions based on the nuclear data for K = 2, with each vertical bar representing an individual.
).For elevation, the resistance cost decreased until approximately 4300 m a.s.l., which might correspond to the relatively flat structure of the Sanetti plateau in the southern region above ~4000 m a.s.l.(FiguresS3c, S5c).The composite resistance surfaces did not show any contribution to genetic differentiation (FiguresS4, S6).As in the partial Mantel tests, we could not identify a contribution of the vegetation index or soil moisture index to genetic differentiation (Table2).4| DISCUSS IONUsing complete mitochondrial genomes and low-coverage nuclear genomes from 77 giant root-rat individuals, we uncovered a clear subdivision of localities into the north and the south of the species range in the Bale Mountains, Ethiopia.Landscape genetic analysis identified topographic barriers (steep slopes and elevation differences) between the two regions as the main drivers of population F I G U R E 3 Analysis of signals of ancient gene flow between individuals WM07 from locality S1 and GG01 from locality S6 and their source group in the south, using D-statistics.Negative D-scores suggest gene flow or recent common ancestry between H1 and H3 relative to H2 and H3, while positive D-scores suggest gene flow or recent common ancestry between H2 and H3 relative to H1 and H3.Statistical significance is indicated by asterisks (*) next to scores, when |Z| greater than 3 (Figure S2), determined by a one-sample Wilcoxon signed rank test.F I G U R E 4 Levels of genetic differentiation between the sampled giant root-rat localities from the north and south of Bale Mountains National Park.Mitochondrial (red, above diagonal) and nuclear (blue, below diagonal) F ST -estimates.* in cells indicates significant differences (p-values < .05),derived by permuting haplotypes between localities for mitochondrial data and by applying a one-sample t-test on F ST -estimates for the nuclear data.subdivision at the landscape scale.Within regions, we did not identify any clear spatial structuring, suggesting a high level of gene flow at the local scale, when topographic barriers are absent.
tochondrial and nuclear genomes was evidenced, for instance, by a high number of substitutions separating the divergent mitochondrial lineages present in each region, or by the nuclear PCA analysis (Figure2, FigureS1).Notably, two giant root-rat individuals sampled in the south were mitochondrially more closely related to their northern counterparts than their source region.This might indicate either recent or ancient landscape-scale gene flow between the regions north and south.However, admixture proportions of nuclear genomes contradict the scenario of recent gene flow, as our findings showed no sign of admixture between regions (Figure2d).Additionally, the D-statistics, specifically testing for ancient gene flow between the two southern individuals and the north, excluded the possibility of ancient gene flow (Figure3).Moreover, these two southern individuals, despite their mitochondrially closer relationship to the northern group, exhibit a large number of substitutions distinguishing them from the northern individuals.Considering these results, in the light of the insights from the nuclear analyses, we propose the clustering of the two southern individuals to the north is due to incomplete lineage sorting(Ballard & Whitlock, 2004;Toews & Brelsford, 2012), a remnant of their shared evolutionary history, rather than recent or ancient gene flow.The phylogeny indicated these two southern individuals as basal to each of two distinct haplogroups found in the north, suggesting the lineages were derived from two distinct divergence events.
soil depth and moisture are likely correlated (deeper soil can store more water,Tromp-van Meerveld & McDonnell, 2006), soil moisture as a proxy for soil availability for giant root-rats may not capture areas of sufficient soil depth.The vegetation and soil moisture indices were derived from Sentinel-2 scenes captured throughout 2021.It is possible that the average values might not fully reveal the specific habitat requirements of the giant root-rat, with regard to its preference for moorlands and wet grasslands with good soil depth.
namics and environmental features that drive the structuring of rangelimited fossorial species.With ongoing environmental changes, it is crucial to utilize this knowledge to safeguard mountain biodiversity and ecosystem functioning.AUTH O R CO NTR I B UTI O N S NF, GM, LO, TW, and DGS designed the research concept and NF, LO, and DGS secured the project funding.VMR, AA, LW, and DGS captured specimens in the field.VMR and AR-I conducted lab work.VMR, AR-I, and MVW analysed the data with contributions from AA, RS, EDL, and DGS for genetic and ecological interpretation.LW generated raster layers for landscape genetic analyses.VMR, EDL, and DGS wrote the manuscript with contributions from all co-authors.All authors gave final approval for publication and agreed to be held accountable for the work carried out within this publication.ACK N O WLE D G E M ENTSThis work was supported by the German Research Council (DFG) in the framework of the joint Ethio-European DFG Research Unit 2358 "The Mountain Exile Hypothesis.How humans benefited from and re-shaped African high-altitude ecosystems during Quaternary climate changes" [FA-925/14-1], [OP-219/10-2] and [SCHA-2085/3-1].
(Groos et al., 2021)ssendorf et al., 2019).Except for this last glacial extent, exposure ages of moraines in the valleys in the northwestern part of the plateau (up to ~100 kya), and stone stripes close to the mountain Tullu Dimtu (up to ~360 kya) could be interpreted in favour of earlier glacial periods(Groos et al., 2021).Possibly, giant root-rat individuals in the south were pushed towards the outer margins of the plateau by glaciers, which increased the separation to the individuals in the north.As the glaciers retreated, colonization of the central plateau from a more southern Late Pleistocene refugium may explain the significantly lower nuclear nucleotide diversity in the central localities (S1 and S2) in comparison with localities in the south east (S3, S5, S6,