DNA barcoding of the mesic adapted striped mouse, Rhabdomys dilectus in the Eastern Cape and KwaZulu- Natal provinces of South Africa

South African small mammals are under-represented in DNA barcoding efforts, particularly from the eastern forested regions of the country. This study reports DNA barcoding of Rhabdomys taxa from previously unsampled parts of the Eastern Cape and KwaZulu-Natal provinces of South Africa. The complete mitochondrial DNA cytochrome oxidase I (COI) gene was sequenced for 101 Rhabdomys sp. individuals from 16 localities from all three main forest groups (coastal, mistbelt, and scarp forests). Molecular data were supplemented with external morphological measurements, including those deemed potential taxonomically diagnostic characters. Findings indicate the area to be inhabited solely by Rhabdomys dilectus chakae. Haplotypes distributed across the three forest groups were separated by shallow sequence divergences ranging from 0.001–0.015 (Kimura 2-parameter model) and displayed very little population genetic structure (FST= 0.071787). Morphological data revealed some regional metric differences in external morphology, but all the head-and-body to tail (HB: tail) ratios match that of R. d. chakae, and consequently, molecular and morphological data are congruent. These data confirm a range extension of R. d. chakae, supporting the utility of COI barcodes in the identification of small mammalian species.


Introduction
The four-striped grass mouse genus Rhabdomys (Thomas, 1916), is distributed throughout southern Africa, occurring in a variety of habitats with ample grass cover (Skinner and Chimimba 2005;Monadjem et al. 2015;Du Toit et al. 2016). Its distribution includes forest edges, where it is a pioneer species following forest clearances (Castley 1997). Rhabdomys is considered a species com-somal data (Rambau et al. 2003;du Toit et al. 2012;Ganem et al. 2020). Of these, R. pumilio and R. dilectus were the first to be recognized based on chromosomal work combined with molecular data (Rambau et al. 2003), with the latter further divided into two subspecies R. d. chakae and R. d. dilectus. Subsequently, further molecular analyses retrieved three paraphyletic lineages of R. pumilio and two monophyletic lineages of R. dilectus, each with distinct geographical distributions (du Toit et al. 2012). Two of these paraphyletic lineages of R. pumilio have been referred to as R. bechuanae and R. intermedius (Ganem et al. 2020). There are contact zones between the R. pumilio, R. bechuanae and R. intermedius that have been identified in the Free State at Sandveld du Toit et al. 2012;Ganem et al. 2020); and at Fort Beaufort in the Eastern Cape (du Toit et al. 2012;Ganem et al. 2020). Other contact zones for R. d. chakae and R. d. dilectus are described from the northern Gauteng Province (Le Grange et al. 2015, Ganem et al. 2020. Altogether, fifteen contact zones for Rhabdomys have been documented in South Africa (Ganem et al. 2020), and while one of these has been identified in the vicinity of the Eastern Cape Province in Fort Beaufort, this has been based on limited sampling despite the species being known to occur there. It remains unclear whether the Eastern Cape population is an extension of R. pumilio, which typically occurs in the southern parts of Eastern Cape Province, or R. dilectus (specifically R. d. chakae), which is associated with most sampled parts of the Eastern Cape and the KwaZulu-Natal provinces.
In the Eastern Cape, there are broad forest group types that extend to the neighboring KwaZulu-Natal and Western Cape provinces: Southern mistbelt, Scarp, and Southern coastal forest groups (Von Maltitz et al. 2003;Cowling et al. 2004;Mucina et al. 2006). South African forests cover a small land surface area (Low and Rebelo 1996), and their distribution is severely fragmented (Geldenhuys 1989;Lawes 1990;Eeley et al. 1999). The south-eastern forests support a high faunal diversity (Minter et al. 2004;Harrison et al. 1997;Skinner and Chimimba 2005;Lawes et al. 2007), that includes a variety of small mammal species that may be classed as either forest dependent or forest associated species (Lawes et al. 2007). Gaining a better understanding of the responses of forest faunal diversity to habitat fragmentation should allow categorization of these species according to these responses (Laurance 1990; Bierregaard and Stouffer 1997;Malcolm 1997;Ranta et al. 1998). In terms of forests, three categories may be identified: interior species, generalist species and edge species . Rhabdomys species are considered generalists and are associated with the forests but are seldom found within forests (du Toit et al. 2012Ganem et al. 2020). Nonetheless, they colonize severely fragmented forest landscapes of eastern South Africa (Ferrante et al. 2017), due to their adaptability to disturbed habitats (sensu Ricklefs and Schluter 1993). As far as we have been able to ascertain, populations of Rhabdomys occurring in and around the forests of the Eastern Cape have not yet been identified using molecular techniques, and therefore it remains uncertain which of the four recognized species occupies these habitats. However, from all the contact zones of the Eastern Cape Province, only three species occur within the province. Rhabdomys pumilio has a continuous distribution from the Eastern Cape Province along the coastline of the Western Cape Province towards the Northern Cape Province. Rhabdomys intermedius is distributed from the Eastern Cape Province hinterland, through the borders of the Western and Northern Cape provinces towards the center of the Northern Cape Province. This leaves the eastern parts of South Africa dominated by R. dilectus, whose range extends towards the north (du Toit et al. 2012;Monadjem et al. 2015;Ganem et al. 2020). Of the two R. dilectus subspecies, R. d. chakae appears to be distributed in the southern parts of KwaZulu-Natal and pockets of the Eastern Cape (Matthee et al. 2018;Ganem et al. 2020), although this requires confirmation by DNA barcoding.
Barcoding provides a rapid way to identify organisms without the need for morphological expertise (Hebert et al. 2003;Hajibabaei et al. 2007), provided that reference material is available. It is possible to assign sampled specimens to species level and also to detect the presence of cryptic species (Meyer and Paulay, 2005;Jacquet et al. 2012). DNA barcoding is a particularly useful method to identify species using the COI short fragment, even for rodents (Borisenko et al. 2008;West 2016;Ramatla et al. 2019), and for cryptic species, e.g., Heller's broad-nosed bat Platyrrhinus helleri Peters, 1866(sensu Clare et al. 2007). One feature of small mammals that makes taxonomic identifications challenging is convergence in external morphology and pelage coloration (Lai et al. 2008;Losos, 2011;Jacobs et al. 2013), making validation using molecular data crucial. Common errors associated with field identifications often result from the absence of referral or comparative specimens; ambiguity due to undescribed species complexes; and verification of the age of the specimen (Borisenko et al. 2008). The benefits of correct taxonomic identification of rodent taxa include recognition of species that are carriers of diseases (Ramatla et al. 2019); an understanding of their role in parasitic transmission cycles and distribution patterns (Engelbrecht et al. 2016); and knowing which rodents are agricultural pests Engelbrecht et al. 2016;Ramatla et al. 2019). To date, several studies have demonstrated the utility of COI in the differentiation of geographic forms of Rhabdomys (du Toit et al. 2012;Grange et al. 2015;Engelbrecht et al. 2016;Coetzer and Grobler 2018;Matthee et al. 2018;Ganem et al. 2020). In the first instance, COI has been used successfully to differentiate museum specimens into two groups -R. dilectus and R. pumilio (Coetzer and Grobler 2018), and more recently R. pumilio group has been divided further to R. pumilio, R. bechuanae and R. intermedius (Ganem et al. 2020). While these studies cover most of the range of the genus in South Africa, there is a noticeable absence of data from parts of the Eastern Cape, particularly from Alexandria to the southern parts of KwaZulu-Natal. Roberts (1951), Coetzee (1970), and Ganem et al. (2020) Smither and Wilson 1979;Rautenbach 1982;Lych 1994;Taylor 1998;Skinner and Chimimba 2005;Ganem et al. 2020). These data suggest that mesic-adapted species exhibit shorter tails (compared to their body size) than their arid-adapted conspecifics. However, the proposed geographically linked morphological variation has yet to be assessed in many Eastern Cape Rhabdomys populations.
Here we present cytochrome oxidase subunit I (COI) barcodes from the Eastern Cape and parts of KwaZulu-Natal forests, to address the distribution gaps highlighted in recent publications (du Toit et al. 2012;Coetzer and Grobler 2018;Ganem et al. 2020). We hypothesize that the COI barcoding would be successful in identifying R. d. chakae as the taxon distributed from Alexandria in the Eastern Cape to Nxumeni state forest in KwaZulu-Natal. Based on the generalist nature of Rhabdomys taxa, we hypothesize additionally that individuals would exhibit genetic relatedness within the sampled distributional range.

Specimen collection and morphological analyses
Animals were collected from three forest types in the Eastern Cape and southern parts of the KwaZulu-Natal Province of South Africa ( Figure 1, Table 1). Small and large folding Sherman traps were placed approximately 2 m apart along traplines in the forest, forest edge, and the grasslands in the surrounding locality. Standard measurements including total length, tail length, hindfoot length (with claw) and ear length, to the nearest 1 mm were taken in the field. Some animals with truncated tails were collected in this study and their measurements were not used; nor were measurements for individuals considered as juveniles or sub-adults. Data were log-transformed and pooled according to forest type and not locality, in order to bolster sample sizes. Descriptive statistics were calculated, and data were screened for normality, skewness, and kurtosis. Owing to uneven sample sizes, and some data sets that did not conform to normal distribution, a Kruskal-Wallis test was used to test for statistical difference between individuals of the three forest types. Statistical procedures were conducted using SPSS 27.0 (IBM Corp. 2020). For molecular analysis, tissue samples were collected from released animals, and voucher specimens, where applicable, were deposited in the Durban Natural Science Museum, South Africa. Supplementary material from the same museum was also used for this study.

Molecular analysis
DNA was extracted using the NucleoSpin (R) Tissue Kit technique (Macherey-Nagel) following the manufacturer's manual. The COI gene fragment was isolated using L1490 (5`-GGT CAA CAA ATC ATA AAG ATA TTG G-3`) and H2198 (5`-TAA CTT CAG GGT GAC CAA AAA ATC A-3`) primers (Folmer et al. 1994). The PCR reactions consisted of 0.8 μM of each primer, 1 unit of Taq DNA Polymerase, 1 X Buffer, 0.2 mM of each dNTP, 2.5 mM of MgCl 2 and distilled water to make up 25 ul reaction volumes. Cycling was performed using a Gene-Amp PCR system (Applied Biosystems). PCR conditions involved an initial denaturing step at 94 °C for 4 min followed by 35 cycles of 30 s denaturation at 94°C, 30 s, annealing at 42.5°C, 30 s extension at 72°C, and then a final extension step at 72°C for 12 min. PCR products were visualized with one percent agarose gel and sequenced using big dye chemistry then sequenced at Central Analytical Facilities of Stellenbosch University.

Sequence data analysis and phylogenetic reconstruction
Sequences generated were edited and assembled manually in BioEdit v7.0.9 (Hall 1999) and aligned using ClustalW (Thompson et al. 2002). Sequences were run on the Basic Local Alignment Search Tool (BLAST, Altschul et al. 1990) at the NCBI website (http://www.ncbi. nlm.nih.gov/blast) for similarity check with published COI sequences of Rhabdomys on GenBank. Sequences were then deposited on the Barcode of Life Data System (BOLD, Ratnasingham and Hebert 2007) through the online interface at www.barcodinglife.org under the Foundational Biodiversity Information Programme (FBIP) in the Foundational Biodiversity Eastern Cape project.
A haplotype dataset was created using DnaSP v6 (Rozas et al. 2017) from aligned sequences and it was subsequently used to generate a network tree using TCS networks (Clement et al. 2000) in PopART (http://popart. otago.ac.nz). Haplotype diversity, population size changes and population demographic expansions (Tajima's D) were also generated using DnaSP v6 (Rozas et al. 2017). Individual sequences were grouped according to the forest type they were collected from, and the p-distances among and within the forest types were calculated using the Kimura 2-parameter model (Kimura 1980); these evolutionary analyses were conducted in MEGA version X (Kumar et al. 2018). To test the significance of population structure, three-way analyses of molecular variance (AMOVA) were done. The genetic differentiation among populations (as defined by forest types) was calculated by pairwise F ST tests. AMOVA and F ST tests were calculated in Arlequin 3.1 (Excoffier et al. 2005) with 1 000 permutations of AMOVA.
All sequences named R. dilectus and R. pumilio retrieved from GenBank were aligned with sequences generated from this study. Phylogenetic trees were reconstructed using Neighbour-joining (NJ), Maximum Parsimony (MP), and Maximum Likelihood (ML) based on the best substitution model which was identified using MODELTEST using MEGA version X (Kumar et al. 2018); support for the nodes was obtained by bootstrap resampling (1000 iterations). GenBank derived COI sequences for Apodemus flavicollis Myomyscus sp., Praomys sp., Rattus rattus, and the sequence of Otomys irroratus (collected from this study) were used as outgroups.
The Barcoding Gap is defined as the difference between the greatest intraspecific variation and the smallest interspecific variation in a group of organisms and is used as the basis for species delimitation (Meyer and Paulay 2005). The mean sequence divergences among specimens were generated using the Kimura 2-parameter (Kimura, 1980) commonly used for barcoding analysis. For estimating the barcoding gap and testing the species hypothesis we used the Automatic Barcoding Gap Discovery (ABGD, Puillandre et al. 2012). Sequences were uploaded on the online interface of ABGD on https:// bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html and analyses were carried out using default settings selecting the Kimura (K80) TS/TV (2.0). Pairwise differences between individual sequences were calculated and their distributions were grouped to form the histogram using the ABGD online interface. In the same way, the number of preliminary species hypotheses (PSHs) obtained for each prior intraspecific divergence were recorded.

Sequence data and BLAST search results
The COI gene fragment (658 bp) was successfully amplified for all specimens (N=101). The nucleotide content of the gene fragment comprised 30.3% of thymine, 25.6% adenine, 27.2% cytosine and 16.9% guanine. Of the 658 bp, only 13.98% was variable (92/658 bp) while 8.5% of sites were parsimony informative and retrieved 56 haplotypes (Table 1). However, when these sequences were combined with all downloaded sequences, the alignment had a few unaligned bases that had to be trimmed. The new alignment contained 280 sequences that were 582 bp long, excluding outgroups. BLAST analysis (Altschul et al. 1990) of our COI sequences attained 99% similarity with R. dilectus, with KC296586.1 top of the list.
Modeltest selected the HKY+G+I (Hasegawa-Kashino-Yano with Gamma and invariant) for the COI dataset. However, because of barcoding analysis preference, phylogenetic reconstructions were undertaken using the Kimura 2-parameter model (Kimura 1980) for NJ and MP. Maximum Likelihood (ML) trees were constructed with both models, and the resulting trees were similar. The sequences downloaded from GenBank include R. bechuanae, R. dilectus (inclusive of R. d. dilectus and R. d. chakae), R. intermedius and R. pumilio. The phylogenetic tree retrieved six clades which were robustly supported with high bootstrap values for ML, MP, and NJ. The me-  Table 1 while the dashed lines indicate the boundaries of each forest type. The relative position of the study area within South Africa is indicated by the rectangle in the insert. sic adapted individuals of R. d. chakae and R. d. dilectus shared a recent common ancestor, while the arid-adapted individuals separated into two paraphyletic clades. One clade contained the individuals of R. pumilio coastal A (individuals ranging from the Eastern Cape Province coast downwards the coastline to Western Cape Province and upwards through the coast of Northern Cape Province in South Africa to Namibian coastal regions) and R. pumilio coastal B (restricted to Eastern Cape around Fort Beaufort). The other clade contained the sister species R. bechuanae (with individuals distributed from Free State Province, South Africa to Namibia) and R. intermedius (individuals from Eastern Cape Province towards the Western Cape Province along with the central parts of these provinces). The insert in Figure 2 indicates all specimens of R. d. chakae.

Rhabdomys dilectus
The mean p-distance among and within forest types for R. dilectus varied from 0.001±0.001 to 0.017±0.003 (Table  3a). The lowest divergence was within the Eastern Cape dune forest grouping (0.001±0.001), while the highest divergence was found amongst Pondoland Coastal scarp forest animals (0.017±0.003). Most observed divergences within the forest types were higher than the divergences Table 1. Summary of sample details including the study localities with corresponding letters (#) matching the map in Figure 1. GPS coordinates, the number of COI (658bp) sequenced individuals and the haplotype numbers are shown, while sample identifiers are shown for voucher specimens and historical museum material only.  (Table 3a). Variation within the population remained high even when the genetic structure was examined between broader forest groups (i.e., Southern coastal forest Southern mistbelt, and Scarp groups).
The median-joining networks for COI (Figure 3) show ed the majority of haplotypes differ by one mutation (i.e., singletons); and those haplotypes were not restricted to forest types or groups. A high number of haplotypes were from the mistbelt and scarp forests, with common haplotypes dominated by individuals from these localities. Of the six common haplotypes, four comprised individuals from both mistbelt and scarp forests, while only two of these haplotypes comprised individuals from coastal forests and either mistbelt or scarp forests. Three of the common haplotypes with individuals from the mistbelt and scarp forests were ancestral haplotypes, with most other haplotypes joined by one Figure 2. A Neighbour-joining (NJ) phylogram illustrating phylogenetic relationships among haplotypes generated from a total of 280 Rhabdomys COI sequences. The dataset includes 101 sequences from this study with other Rhabdomys sequences downloaded from GenBank. The included outgroups were two sequences of Praomys species and Myomyscus species, and a single sequence of Rattus rattus, Apodemus flavicollis and Otomys irroratus. The tree is condensed, while the insert shows R. d. chakae which includes all samples from this study. The Neighbour-joining Maximum Likelihood (ML), and Maximum Parsimony (MP) analyses were drawn using the Kimura 2-Parameter model (Kimura 1980). Bootstrap resampling support (1000 iterations) are shown for uncut nodes, respectively. Evolutionary analyses were conducted in MEGA version X (Kumar et al. 2018). mutation from these haplotypes. Although there were no clear patterns of forest segregation indicated by the haplotypes, there were few haplotypes that are individually divergent from the others. A total of 56 haplotypes (from n=101) resulted in a high haplotype diversity of 0.972 (Table 1).

Barcoding Gap among Rhabdomys species
The groupings as retrieved above by the phylogenetic tree were used for calculating divergences between and within Rhabdomys taxa (Table 3b). The results revealed that the nucleotide difference was always smaller within all the groups (between 0.004±0.001 and 0.023±0.004) than the lowest nucleotide difference between groups (which is between R. pumilio coastal A and coastal B groups 0.050±0.002). The highest nucleotide difference within a group was for R. d. dilectus (0.023±0.004); this nucleotide difference was still lower than the lowest nucleotide difference between groups: between R. pumilio coastal A and coastal B groups (0.050±0.002). There was a notable lower nucleotide difference within R. d. chakae (0.008±0.002), compared to the nucleotide difference between R. d. chakae and the other groups. The low nucleotide differences between groups were between R. pumilio coastal A with R. pumilio coastal B (0.050±0.002); between R. d. chakae and R. d. dilectus (0.053±0.001); and between R. bechuanae and R. intermedius (0.079±0.001). The ABGD analyses (Figure 4) Figure 1). The number of mutations separating haplotypes is indicated using hatch marks. The network was drawn in PopArt.   Table 2. Morphological analyses of 78 Rhabdomys dilectus specimens across different forest types. Standard measurements including total length (TL), tail length (T), hindfoot length (with the claw, HF) and ear length (E), to the nearest 1 mm were taken in the field; the head and body length (HB) was calculated as the difference of tail and the total length of the animal. Sample sizes are presented in parentheses. Descriptive statistics (mean, standard deviation) and Kruskal-Wallis tests were conducted in SPSS 27.0 (IBM Corp., 2020).  (Rambau et al. 2003;du Toit et al. 2012;Ganem et al. 2020). With the exception of R. d. dilectus (with a variation slightly higher 0.023±0.004), there was very low genetic variation within populations of each of the five groups (between 0.005±0.002 and 0.009±0.002), suggesting that R. d. chakae could be an independent species. However, the low variation between R. pumilio coastal A and R. pumilio coastal B (0.050±0.002); and between R. d. chakae and R. d. dilectus (0.053±0.001), which comprise the lowest variation between groups, suggest that these groups are subspecies or lineages (R. dilectus and R. pumilio) respectively. This is supported by the finding that divergences for small mammals with more than one distinct lineage can have more than 3 % intraspecific variation (only up to a maximum of 5.4 %), and only a few species without known lineages show more than 2 % intraspecific variation (Borisenko et al. 2008). To support this the ABGD (Puillandre et al. 2012), preliminary species hypothesis (PSH) suggested that six distinct groups should be considered as species, and the distance phylogenetic tree revealed six groups and each group had at least 99% bootstrap support.

Mistbelt
Our combined (sequence divergence, ABGD and phylogenetic analyses) findings are largely congruent, and suggest six groups may be assigned to two sub-species of R. dilectus; two lineages of R. pumilio; and two species R. bechuanae and R. intermedius. These findings are consistent with the previous findings relating to this genus (du Toit et al. 2012;Ganem et al. 2020). But our study provides much-needed barcoding data for Rhabdomys from the Eastern Cape Province which was not provided by previous studies (Castiglia et al. 2011;Meynard et al. 2012;Ganem et al. 2012;Ganem et al. 2020). Our study documents the distribution of this taxon in the unsampled parts of the Eastern Cape and southern KwaZulu-Natal Table 3a. The mean p-distances among and within (in bold) forest types in the Eastern Cape and KwaZulu-Natal provinces of South Africa for COI (658bp) for Rhabdomys dilectus. The number of base substitutions per site, derived from the mean of all sequence pairs between groups, is shown with the standard error estimate(s) next to each. Analyses were conducted using the Kimura 2-parameter model (Kimura 1980). Evolutionary analyses were conducted in MEGA version X (Kumar et al. 2018 Table 3b. The number of base substitutions per site in Rhabdomys taxa. These were derived from the mean of all sequence pairs between and within species, with standard error estimate(s) shown next to each. Analyses were conducted using the Kimura 2-parameter model (Kimura 1980). The analysis involved 280 nucleotides in the COI dataset which had a total of 582 bp. Evolutionary analyses were conducted in MEGA version X (Kumar et al. 2018 Morphology Ganem et al. (2020) suggested that morphological comparison within Rhabdomys is unreliable, because whenever there is more than one collector, their measuring criteria may be different. But despite this reservation, the head-and-body to tail ratio has been used to separate the mesic group from the arid group (Ganem et al. 2020). The analysis involved controlling for mass and data collector effect and found that tails of the mesic species R. dilectus (R. d. dilectus and R. d. chakae) were, in general, shorter than those of their arid-adapted counterparts. Our study revealed possible clinal variation in R. d. chakae external morphology, with individuals from more temperate climates (mistbelt) slightly larger or longer than their coastal conspecifics. Animals recorded from coastal forests appeared to be the shortest of the three groups investigated. However, we found no significant difference in the ratio of the head-and-body to tail across the forest types which averaged 1.11±0.03 mm, matching the ratio of R. d. chakae from Ganem et al. (2020). This confirms that all specimens in this study are indeed R. d. chakae, supporting the taxonomic assignments based on DNA barcoding in this group.

Phylogeography
The reticulate tree (Figure 3) indicates that this species has haplotypes comprising individuals from all over the study area. It should be noted that most of the shared haplotypes include individuals from the mistbelt and scarp forests, including almost all single haplotypes. The distribution and assemblage of forest animals in South Africa have been heavily influenced by climatic fluctuations and species recolonization following these fluctuations. However, forest-associated species such as R. d. chakae were less affected by such events, as they were able to survive in adjacent habitats (Lawes et al. 2007). The high number of single haplotypes and high haplotype diversity within R. d. chakae may be attributed to very high reproductive rates: an average of two litters per season producing up to 11 offspring per litter (Taylor 1998;do Toit et al. 2016). There was no specific structure detected with the AMO-VAs, and the phylogenetic tree (insert in Figure 2) has robust bootstrap support. This indicates that our sampled individuals are indeed one species, with no indications of lineages in this species, and exhibit a low variation of less than 1% within this group (0.008±0.002, its generalist nature as it utilizes a wide range of resources (Brown 1984). Although a wide range of rodents are good candidate taxa for phylogeographical analysis (Bryja et al. 2014), it is important to note that rodents have different habitat responses and sensitivity to environmental change . Generalist species, such as the subspecies R. d. chakae are hypothesized to conform to the predictions of meta-population theory rather than island biogeography (MacNally and Bennet 1998). Indeed, the genetic analysis of R. d. chakae supports this postulation, as no phylogeographic structure was detected. These data indicate that there is connectivity between the populations in different forests and reflects the ability of the species to survive in disturbed, semi-urban and agricultural landscapes ). The distribution of R. d. chakae may have been limited by the Bedford gap (Lawes 1990;Lawes et al. 2000;Lawes et al. 2007), which has been shown to limit the distribution of other small mammal species (Willows-Munro 2008, Abrahams 2018. There are no records of the subspecies beyond this gap: south-west of the Bedford area there is a population of R. pumilio extending to the Western Cape Province (du Toit et al. 2012). The recorded contact zones in the Eastern Cape Province are in the Bedford area (du Toit et al. 2012;Ganem et al. 2020) and may be attributed to recent colonization (as Tajima's D indicated) that the population of R. dilectus from the sampled area has undergone recent expansion.

Conclusions
The utility of DNA barcoding in identifying Rhabdomys taxa has been supported. This study focused on supplying information on the distribution of the taxa in the Eastern Cape and southern parts of KwaZulu-Natal provinces and confirmed that all samples included in this study were of R. d. chakae. These findings are consistent with previous revisions of this genus (du Toit et al. 2012, Ganem et al. 2020, suggesting four species with two subspecies of R. dilectus and two geographic and genetic distinct lineages of R. pumilio. Our molecular and morphological data support the delineation of R. d. chakae and R. d. dilectus . There may be intraspecific morphological clinal variation within the populations studied here, which merits further investigation. Our findings support the notion that although the ecological requirements of species may differ geographically, generalist species tend to show connectivity between subpopulations within fragmented habitats.