Ethical clearance was obtained from the Interfaculty Animal Ethics Committee at the University of the Free State, South Africa (Ethics number: UFS-AED2018/0075).

A farmer from the Cradock district in the Eastern Cape reported the sighting of a T. swinderianus carcass on his property. A tail clipping of the carcass was submitted to the Genetics laboratories at the University of the Free State for future analysis. This specimen will be referred to as CR_CDK_01_EC. The sample was stored in 96% ethanol at –20°C. This species has never been seen in this region. The farm is situated in a valley nestled within the Winterberg Mountain range. The habitat consists of mixed vegetation with the Dry Highveld Grassland Bioregion on the mountain plateau and the Sub-Escarpment Grassland Bioregion (Mucina and Rutherford 2006) in the valley with a seasonal river draining into the Tarkariver system. The banks of these river systems consist of a mixture of tree and shrub species, including Acacia (Vachellia) karroo, Searsia sp. and Lycium sp. The Tarkariver system in turn runs into the Great-Fish River, which enters the Indian Ocean east of Grahamstown.

Additional samples were sourced from South African museum collections to use in the population genetic and phylogenetic analyses (n = 30; Table S1). The museum specimens were sourced from four museums, namely the National Museum (Bloemfontein), Ditsong National Museum of Natural History (Pretoria), Durban Natural Science Museum (Durban) and Amathole Museum (King William’s Town). The museum samples consisted of either dried skins or skull scrapings. Further details on each sample can be obtained in Table S1. Successfully sequenced specimens were grouped according to the provincial region of origin (Fig. 1).

All DNA extractions were performed with the Purelink Genomic DNA Kit (Life Technologies, Carlsbad, CA). The manufacturers protocol was followed for DNA extraction from animal tissues, with additional steps to ensure sufficient DNA quantity and quality as reported by Coetzer and Grobler (2019). A NanoDrop Spectrophotometer ND-1000 was used for DNA quantification. A dedicated ancient/museum DNA laboratory was not available for DNA extraction. All work surfaces and equipment were decontaminated using a 1.25% hypochlorite solution and then a 70% ethanol solution before and after each set of DNA extractions.

Partial segments of two mitochondrial regions were targeted for this study. Amplification of a ~650 bp fragment of the cytochrome c oxidase I (COI) gene was performed using primers designed by Gaubert et al. (2014). A fragment of ~500 bp was amplified for the mitochondrial D-loop region using primers reported by Adenyo et al. (2013). PCR setup was performed in a DNA free laboratory, and all working surfaces were decontaminated using a 1.25% hypochlorite solution, prior to setup. Each 10 µl PCR reaction consisted of: ~50 ng template DNA, 5 µl Ampliqon TEMPase Hot Start 2 × Master Mix, 0.4 µM of each primer, 0.5 mg/ml BSA and PCR grade dH₂O to make up the 10 µl. The PCR cycling conditions included an initial denaturation step at 95°C for 3 min, followed by 35 cycles of 95°C for 30 s, 58°C for 40 s and 72°C for 30 s, with a final extension step of 72°C for 1 min. Negative PCR controls were included with each PCR reaction batch to monitor the occurrence of possible contamination.

Amplification success was assessed on a 1% agarose gel. The ExoSAP-IT PCR Product Clean-up kit (Affymetrix Inc., CA, USA) was used for all PCR clean-up procedures. The ABI PRISM BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems Division, Perkin-Elmer, Foster, CA) was used for all sequencing PCR reactions. Sequencing PCR clean-up was performed with the BigDye XTerminator Purification Kit (Applied Biosystems, CA, USA), followed by sequencing analysis on an ABI 3500 Genetic Analyzer.

Previously reported COI (n = 23; accession numbers: KJ192912-KJ192913, Equatorial Guinea; KJ192933– KJ192949, Ghana; KJ192950–KJ192953 and KJ192955, Nigeria) and D-loop (n = 26; accession numbers: AB675385–AB675410, Ghana) sequences were downloaded from GenBank for downstream analyses (Adenyo et al. 2013; Gaubert et al. 2014). Sequence assembly and editing of sequences obtained from the current study were performed with GENEIOUS R9 software (Kearse et al. 2012). All sequences generated from this study were deposited in GenBank (accession numbers: COI, OP121188–OP121208; D-loop, OP121209–OP121231; Table S1). The sequence alignments were performed using the CLUSTALW (Thompson et al. 1994) add-in for GENEIOUS. Summary statistics (number of haplotypes (h), haplotype diversity (Hd), number of polymorphic sites (S), nucleotide diversity (π), number of synonymous sites and number of non-synonymous sites) were estimated in DNASP v5.10.01 software (Librado and Rozas 2009).

The most optimal nucleotide substitution model for haplotype alignments of COI and D-Loop were identified prior to phylogenetic tree construction, using the Akaike information criterion (AIC; Akaike 1974) as implemented in JMODELTEST v.2.1 (Darriba et al. 2012). These were identified as HKY +G for both COI and D-loop. Phylogenetic analyses for each gene region were performed via Bayesian inference (BI) in MRBAYES v3.2 (Ronquist et al. 2012). Fukomys damarensis and Hystrix indica COI and D-loop sequences were used as outgroups (Table 1). Ten million generations, a sampling frequency of 1 000 and a burn-in of 0.25 was selected for the BI analysis of both genes. The average standard deviation of split frequencies was used to assess run convergence (split frequencies < 0.01). A minimum spanning network was constructed for each gene region using the program PopArt (Leigh and Bryant 2015). Between group mean pairwise Kimura 2-parameter model (K2P; Kimura 1980) genetic distance measures were calculated for each gene region using MEGA 11 (Tamura et al. 2021). This analysis was performed to assess the level of genetic distance between the South African and West African T. swinderianus groups. For comparative purposes, an additional 11 Hystricomorpha and 3 Sciuromorpha species were included (Table 1). The K2P analysis were run for 1 000 bootstraps.

A molecular clock analysis was performed in BEAST2 (Bouckaert et al. 2014) to improve our understanding of the phylogenetic history of T. swinderianus using the available data. A partitioned dataset was created, using the COI and D-loop sequences from the South African specimens. Additional COI (Gaubert et al. 2014) and D-loop (Adenyo et al. 2013) sequences were downloaded from GenBank according to geographical region of origin. The outgroup taxa comprised of 14 species, representing four Hystricomorpha families and one Sciuromorpha family as outgroups (See Table 1 for details). The most optimal nucleotide substitution model for the COI data set was identified as GTR +I +G and for the D-loop datasets as HKY+G using JMODELTEST. The lognormal relaxed-clock model, with the calibrated Yule process tree prior was selected. Six calibration dates were used from several sources, including the crown Thryonomyidae group estimate from Kraatz et al. (2013) at 5.6 to 2.6 Mya (Table 2). The simulation was run for 20 million generations, with a tree sampling frequency of 20 000. Simulation stationarity was assessed using the effective sample size (ESS) values obtained from TRACER v.1.7.1 (Rambaut et al. 2018). The maximum clade credibility (MCC) tree was estimated using TREEANNOTATOR v.2 (Bouckaert et al. 2014) after discarding 10% of the sampled trees. The MCC tree was assessed and divergence dates were retrieved using FIGTREE v.1.4.0 (Rambaut 2014).

Due to the small sample size of this preliminary study, all genetic diversity values should be taken with caution. The samples from the KwaZulu-Natal group showed the highest genetic diversity values for both loci (COI, Hd = 0.844, π = 0.003; D-loop, Hd = 0.972, π = 0.006) compared to the other South African regions. Assessing the South African sample set as a whole, it was observed that the COI haplotype diversity (Hd = 0.724) is only slightly lower than that observed for Ghana (Hd = 0.814; Table 3). Six haplotypes were observed among the 21 specimens, with all six haplotypes occurring within the KwaZulu-Natal group. The overall South African D-loop genetic diversity estimates were much lower than observed for Ghana, although there is a large sample size difference to consider. Nine D-loop haplotypes were observed for the South African group, with eight haplotypes observed in KwaZulu-Natal.

The Bayesian inference (BI) trees for COI and D-loop both identified two clearly separated clades, dividing West African and South African specimens (posterior probability = 1; Fig. 2). The pairwise K2P analyses showed higher mean genetic distances between the South African and West African T. swinderianus groups (COI K2P = 0.112, SE = 0.015; DL K2P = 0.117, SE = 0.017), than what was observed for other species pairs tested. All three pairwise genetic distances calculated for the representative Hystrix species at both COI and DL were substantially lower that what was observed for the two T. swinderianus groups. A similar result was also observed for the more distantly related Marmota sp. included in this analysis (Table S2).

Figure 2. 

The phylogenetic trees reconstructed by Bayesian inference for (a) the COI and (b) the D-loop regions. Sequences from South African formed a distinctively unique clade separate from the western African sequences. Posterior probability values for the main clades are provided next to each node. The haplotypes representative of the extralimital Eastern Cape specimen are shown in bold type.

No clear geographic structuring within South Africa was observed within either sequence datasets. From the haplotype network analysis, it can however be observed that there are four COI haplotypes unique to the KwaZulu-Natal group (Fig. 3a). One COI haplotype was shared between the KwaZulu-Natal and Eastern Cape groups, and another was shared by all three groups. A similar trend was observed for the D-loop dataset, with six unique haplotypes recorded for KwaZulu-Natal. One unique haplotype was identified for the Free State group, containing 17 mutation. No COI sequence was obtained for this specimen. One D-loop haplotype was present in all South African groups, with another shared between the KwaZulu-Natal and Eastern Cape groups (Fig. 3b).

Figure 3. 

The minimum spanning networks constructed from (a) COI and (b) D-loop sequences generated from the current study and source from GenBank (see text for details). Distinct southern and western African clades can be observed for both mitochondrial loci.

Specimen CR_CDK_01_EC collected from the Eastern Cape, neatly grouped with the rest of the South African T. swinderianus specimens confirming the species identity. The haplotype network analyses, placed specimen CR_CDK_01_EC within the main haplotype cluster for both loci (Fig. 3). The COI cluster contained all the Free State specimens and one KwaZulu-Natal specimen. The D-loop cluster containing specimen CR_CDK_01_EC consists of specimens from all regions sampled, with all but one Free State specimen also in this cluster. No clear phylogenetic position within the South African clade was observed during the BI analyses.

Divergence date estimates are provided in Table 4, with the letters on Fig. 4 corresponding to the node names in the table. The most recent common ancestor (MRCA) of Phiomorpha rodents were estimated at 33.935 million years ago (Mya) (95% highest posterior density (HPD) = 29–40.991 Mya). The MRCA of Thryonomys was calculated at 2.919 Mya (95% HPD = 2.6–4.154 Mya), which was only older than the date estimated for the split between Hystrix africaeaustralis and H. cristata (0.308 Mya; 95% HPD = 0.047–0.653) (Table 4). The MRCA of the South African T. swinderianus clade was estimated at 0.199 Mya (95% HPD = 0.069–0.373 Mya) and the MRCA for the West African clade was calculated as slightly older at 0.277 Mya (95% HPD = 0.122–0.476 Mya).

Figure 4. 

The maximum clade credibility (MCC) tree obtained from BEAST from a concatenated COI and D-loop dataset successfully differentiated between the Rodent suborders used in the current study. The nodes are labelled according to the node names in Table 4. Myo = Myomorpha; Sci = Sciuromorpha; Lag = Lagomorpha

The placement of specimen CR_CDK_01_EC with specimens from the Free State group was unexpected, as the Grahamstown region which is the excepted southern boundary of the T. swinderianus distribution is but 120 km south-east from this site. It was thus originally thought that cane rats travelled along the Great Fish river system, which enters the Indian ocean east of Grahamstown, and travelled to where the specimen was found in the Cradock district (Figure 5). The haplotype network, however, point to a possible migration event southward from the Free State Province. This is feasible since there have been an anecdotal report from van der Merwe and Avenant (2004) of a T. swinderianus sighting along the Leeu River valley, which is an tributary of the Caledon river in the eastern Free State. The Caledon river in turn flows into the Orange river. This region is known for large agricultural developments utilising irrigation from the Gariep dam situated in the Orange River (Nieuwoudt et al. 2004). Further south of the Gariep dam are several crop fields fed from either borehole irrigation or water from the Orange-Fish River irrigation system (W.G. Coetzer, personal observation). It can therefore be argued that cane rats could have migrated south from the eastern Free State along these river systems and patches of cropland, eventually ending up along the Fish River system near Cradock.

Figure 5. 

An elevational map depicting the Great Fish River (pink) system in relation to the new extralimital sample from the Eastern Cape Province. The Orange River system is shown in orange, with the location of the Gariep Dam also indicated in blue.

The use of T. swinderianus as a human protein source, from either bushmeat or captive-bred sources, could also potentially contribute to the occurrence of these animals outside their known distribution ranges. It is suggested by Avenant et al. (2016) that the observation of this species in the central Free State Province could be linked to human mediated translocations.

It is clear that T. swinderianus populations are continuing to expand their range in a westerly direction in South Africa, corroborating a previous report by van der Merwe and Avenant (2004). van der Merwe and Avenant (2004) provided the first reports of this westerly range expansion of T. swinderianus in South Africa, warning of the possible negative effects this could have on agricultural practices in the region. This can further be observed from museum collections outside of the current distribution range in KwaZulu-Natal (Figure 1). Following these authors’ recommendations, it is important for conservation management authorities to closely monitor this range expansion. These rodents are known to cause damage to crops across Africa (Ewer 1969; van der Merwe and Avenant 2004; van Zyl and Delport 2010), therefore appropriate human-animal conflict solutions should be developed before the problem become too big. Additional sampling at a fine scale across the distribution range could provide better insights into the phylogeographic structure of this species, shedding more light on possible range expansion and colonisation events. The addition of nuclear as well as adaptively linked loci would further benefit future studies.

This study is the first to report on T. swinderianus genetic diversity and phylogenetic structure from South Africa. The high level of genetic diversity, and the occurrence of unique haplotypes, in the eastern regions of South Africa supports an ancient colonisation event from East Africa southwards into South Africa. Similar trends have been observed in other mammal species in southern Africa including African mole-rats (Heliophobius sp.; Faulkes et al. 2004), samango monkeys (Cercopithecus sp.; Lawes 1990; Linden et al. 2020), vervet monkeys (Chlorocebus spp.; Turner et al. 2016), and greater kudu (Tragelaphus strepsiceros; Jacobs et al. 2022) to name a few. The small sample size of the current study could have attributed to the lack of significant genetic structuring within the South African distribution range.

The haplotype diversity levels observed for the South African group was only slightly lower than that observed from data obtained from Adenyo et al. (2013) and Gaubert et al. (2014) for west African populations. The absence of any shared haplotypes for either the COI or D-loop regions between western and southern Africa, could point to an ancient common ancestry and independent evolutionary processes driving mutation rates in western vs southern Africa. This was further supported by the high K2P value observed between the two T. swinderianus lineages. This observation points to possible separate colonization events into western and southern Africa from eastern or north-eastern Africa, as fossil records from east Africa suggest this region being the origin for the genus (Winkler 1992; Winkler et al. 2010; Kraatz et al. 2013).

The topology of the maximum clade credibility (MCC) tree obtained from BEAST reflects the known higher taxonomic groupings of the Rodentia taxa used, as compared to other reports (Fig. 4; Upham and Patterson 2015; Swanson et al. 2019). It is worth noting that several contradicting phylogenetic topologies for the Hystricomorpha families have been reported in the past, which for one could be an artefact of the molecular marker used (Huchon et al. 2000; Honeycutt et al. 2007; Heritage et al. 2016; Swanson et al. 2019).

The divergence dates estimated here are, however, in line with previous studies. The accuracy of the estimated Thyronomys divergence dates were assessed by comparing estimates of well-known taxa to the literature. The MRCA date estimated for Phiomorpha in this study is slightly younger than estimates reported by Upham and Patterson (2015), although these estimates do overlap with the confidence intervals calculated in that study. The Phiomorpha date calculated in the current study further overlaps with the confidence intervals of estimates from Sallam et al. (2009), Upham and Patterson (2012) and Huchon et al. (2007). The divergence estimates at lower taxonomic ranks were generally much younger than previously published. An example is the divergence date for the Atherurus/Hystrix split (4.38; 95% HDP = 3.582-5.656). This date estimate is relatively younger than what has been reported by others (8.6-13.6 Mya; Sobrero et al. 2014; Upham and Patterson 2015; Kumar et al. 2017). These lower divergence estimates could be an artefact of the markers used, and should be taken into consideration when viewing the results (Mueller 2006; Swanson et al. 2019). Additional markers could provide more reliable results.

The divergence time estimates for the T. swinderianus clades from the current study is younger than the divergence date estimations by Kraatz et al. (2013), based on thryonomyid fossil records from the United Arab Emirates and assessments of previous publications. Early Thryonomys sp. fossils identify the origin of the genus in east Africa from the late Miocene (6 – 5.6 Mya) (Haile-Selassie et al. 2004; Manthi 2007; Wesselman et al. 2009). The estimated divergence date for the South African T. swinderianus clade indicate the occurrence of a common ancestor during the late Pleistocene. The occurrence of Thryonomys sp. fossils in southern Africa were observed from the late Pleistocene at 0.126 to 0.0117 Mya (Brain 1981; Winkler et al. 2010; ICS 2022). This is slightly younger than the molecular date estimated from this paper, although the 95% HPD does overlap with the fossil ages. The mid- to late Pleistocene consisted of oscillations of wetter and drier conditions in Africa, leading to possible habitat fragmentation and the formation of unique habitats (Steele 2013; Hoag and Svenning 2017). Eastern Africa, specifically, underwent severe drought conditions starting in the middle Pleistocene (Hoag and Svenning 2017). This in turn would lead to population isolation and differentiation due to reduce gene flow and genetic drift. Further studies on the phylogenetic structure of T. swinderianus utilising specimens from across Africa would provide further details on the evolution of Africa’s second largest rodent. Such information could be useful for conservation authorities to establish conservation plans in areas where the species is experiencing anthropogenic pressures which could have long-term detrimental effects.