Relic populations of Fukomys mole-rats in Tanzania: description of two new species F. livingstoni sp. nov. and F. hanangensis sp. nov.

Previous studies of African mole-rats of the genera Heliophobius and Fukomys (Bathyergidae) in the regions of East and south central Africa have revealed a diversity of species and vicariant populations, with patterns of distribution having been influenced by the geological process of rifting and changing patterns of drainage of major river systems. This has resulted in most of the extant members of the genus Fukomys being distributed west of the main Rift Valley. However, a small number of isolated populations are known to occur east of the African Rift Valley in Tanzania, where Heliophobius is the most common bathyergid rodent. We conducted morphological, craniometric and phylogenetic analysis of mitochondrial cytochrome b (cyt b) sequences of two allopatric populations of Tanzanian mole-rats (genus Fukomys) at Ujiji and around Mount Hanang, in comparison with both geographically adjacent and more distant populations of Fukomys. Our results reveal two distinct evolutionary lineages, forming clades that constitute previously unnamed species. Here, we formally describe and designate these new species F. livingstoni and F. hanangensis respectively. Molecular clock-based estimates of divergence times, together with maximum likelihood inference of biogeographic range evolution, offers strong support for the hypothesis that vicariance in the Western Rift Valley and the drainage patterns of major river systems has subdivided populations of mole-rats. More recent climatic changes and tectonic activity in the “Mbeya triple junction” and Rungwe volcanic province between Lakes Rukwa and Nyasa have played a role in further isolation of these extra-limital populations of Fukomys in Tanzania.


INTRODUCTION
African mole-rats of the family Bathyergidae are subterranean rodents that occur throughout sub-Saharan Africa (Faulkes & Bennett, 2013), with much of their range subdivided by the Great Rift Valley. They have been widely studied as a result of the variation in their social and reproductive strategies, and comparative studies have become crucial in this respect (Allard & Honeycutt, 1992;Faulkes et al., 1997;Faulkes & Bennett, 2013).
More recently, the naked mole-rat (Heterocephalus glaber) has also emerged as a model species for the study of longevity and cancer resistance (Gorbunova et al., 2014). Hence a clear understanding of their taxonomy, biodiversity and phylogenetic relationships has become increasingly important, not least because they are a speciose group, but also because there are a number of genetically unique, disjunct populations that are limited in their distributional range (Faulkes et al., 2004;Ingram, Burda & Honeycutt, 2004;Van Daele et al., 2007a;Van Daele et al., 2007b). Historically, systematics of the group has been problematic, because cryptic diversity is prevalent in mole-rats due to convergent morphological evolution of a phenotype adapted to the subterranean niche. However, molecular phylogenies based on both nuclear and mitochondrial genes have produced congruent trees (e.g., Allard & Honeycutt, 1992;Faulkes et al., 1997;Walton, Nedbal & Honeycutt, 2000;Huchon & Douzery, 2001;Faulkes et al., 2004;Ingram, Burda & Honeycutt, 2004;Van Daele et al., 2007a;Van Daele et al., 2007b).
Plate tectonics and the formation of the Great Rift Valley have played a central role in the adaptive radiation and distribution of the Bathyergidae, particularly among mole-rats of the genera Heliophobius and Fukomys (Faulkes et al., 2004;Faulkes et al., 2010;Faulkes et al., 2011). These taxa are distributed widely, with virtually all members of the genus Fukomys occurring in locations west of the Western (Albertine) and Southern Rift Valleys from northern South Africa, through south-central Africa to Uganda and Sudan. South-central Zambia in particular is a hot-spot for species/karyotypic diversity in Fukomys, possibly as a result of changing patterns of drainage over geological time (Van Daele et al., 2004;Van Daele et al., 2007a;Van Daele et al., 2007b). Disjunct populations are found in Ghana, Cameroon and Nigeria, and a small number of isolated populations occur east of the Rift Valley in Tanzania, where the silvery mole-rat Heliophobius is widespread and the predominant bathyergid rodent (Faulkes et al., 2011). Faulkes et al. (2010) investigated a number of populations of Fukomys (or Cryptomys sensu lato) in Tanzania in an attempt to clarify their taxonomic status and to confirm the nomenclature proposed by the earliest reports published by Allen & Loveridge (1933). The latter originally described a new taxon (Cryptomys hottentotus occlusus) from Kigogo in south-western Tanzania, interpreting it as a locally adapted form of Cryptomys hottentotus whytei (Fukomys whytei stricto sensu; Van Daele et al., 2007a), which is geographically the closest in distribution to C. h. occlusus. Allen and Loveridge also report catching F. whytei (stricto sensu) from further north at Ujiji. An additional two, more distant locations (Mount Hanang and Liwale), were later recorded for C. h. occlusus by Swynnerton & Hayman (1951) in their checklist of Tanzanian mammals. The study by Faulkes et al. (2010) concluded that Fukomys whytei constitutes a clear phylogenetic species, supporting the ''whytei'' clade described by Van Daele et al. (2007a), and that C. hottentotus occlusus (sensu Allen & Loveridge, 1933) should be subsumed into F. whytei or, at most, considered a subspecies. With regard to animals sampled from populations at Liwale and Hanang, the former were found to be Heliophobius rather than Fukomys (Faulkes et al., 2011), while genetic analysis of two mole-rats from Hanang appeared to constitute a previously unrecognised species (Faulkes et al., 2010). At the time it was not possible to obtain samples from the remaining sites at Ujiji described by Allen and Loveridge.

Analysis of mitochondrial DNA sequences
Sequences were aligned manually for analysis using Mesquite version 3.03 (Maddison & Maddison, 2014) and phylogenetic relationships investigated using parsimony, maximum likelihood and Bayesian approaches. These methodologies build evolutionary trees using different underlying assumptions and algorithms, and it is usual to investigate the data with all three approaches. Parsimony aims to produce a tree that reflects the minimum number of evolutionary changes needed to explain the data, while likelihood and Bayesian methods take a probabilistic approach based on models of molecular evolution. Maximum Likelihood fits of 24 different nucleotide substitution models were used to establish the evolutionary model most appropriate for the data (from Hierarchical Likelihood Ratio  tests), and these parameters were then used in subsequent analyses, where appropriate. Maximum likelihood and parsimony analyses were undertaken and phylogenetic trees and genetic distances among haplotypes based on nucleotide sequences constructed using MEGA 6 (Tamura et al., 2013). Maximum likelihood was conducted using the heuristic search option, with initial tree(s) for the search obtained automatically by applying the Maximum Parsimony method. For maximum parsimony we used the min-mini heuristic algorithm with a search factor of 1 with gaps treated as missing data and eliminated from the analysis. Bootstrap analysis was conducted with 100 replicates of the dataset. Bayesian phylogenetic analysis was undertaken using BEAST v1.8.2 (Drummond & Rambaut, 2007;Drummond et al., 2012). Following the molecular clock likelihood ratio test performed using MEGA 6 (Tamura et al., 2013) to establish the correct molecular clock model, the null hypothesis of equal evolutionary rate throughout the tree was rejected (likelihood ratio = 9.92; P > 0.001). Thus an uncorrelated relaxed molecular clock model (Drummond et al., 2006) and a Yule tree prior (the most suitable for interspecies comparisons) were selected in BEAST, and an HKY model of molecular evolution. The molecular clock rate was calibrated by assuming a divergence time of 10-11 Mya for the common ancestor of Cryptomys/Fukomys (the ingroup in this study), and these divergence times for the ingroup were input as a prior with upper (11) and lower (10) limits. This calibration has been previously used by Ingram, Burda & Honeycutt (2004), Van Daele et al. (2007a) and Faulkes et al. (2010), and was based on a timing of 19 Mya for the divergence of the Heliophobius lineage within the bathyergid family tree, and the occurrence of the fossil Proheliophobius (Lavocat, 1973). After initial data exploration with independent chains we implemented a final run having a chain length of 30,000,000, sampling output every 30,000 iterations. The first 300,000 trees (10%) were discarded during burn-in. Mixing and convergence of MCMC chains generated by BEAST were investigated and checked using Tracer v1.6.0 to ensure sufficient iterations and sampling were performed before samples from the posterior distribution of trees were summarized using Treeannotator v1.8.2, and trees drawn using FigTree v.1.4.2 (Drummond et al., 2012).
Each distinct haplotype obtained from the two geographical locations (Ujiji and Hanang) were included in all phylogenetic analyses, together with the published sequences representative of the main clades of Fukomys (Van Daele et al., 2007a;Faulkes et al., 2010), and two other bathyergid mole-rats as outgroups: Cryptomys hottentotus hottentotus and Heliophobius emini (Faulkes et al., 2011). In addition, a previously unpublished Fukomys sequence from Ghana was included as another example of an extralimital, but geographically distant population.

Phylogeographic analysis
Phylogeographic analysis and hypothesis testing of the ingroup clade Cryptomys and Fukomys was performed using likelihood analysis of geographic range evolution, utilizing the dispersal-extinction-cladogenesis (DEC) model implemented in the program Lagrange v.2.0.1 (Ree & Smith, 2008). Within the Lagrange DEC model, the cladogenesis parameter remains fixed, while the other two parameters (dispersal/range expansion and extinction/range contraction) are estimated. Seven geographical areas were defined, based on the distribution of extant clades: south of the Zambezi River (SZ); south of the Limpopo River (SL); west of the Kafue River (WK); east of the Kafue River (EK); Mbeya Triple Junction (MTJ); east of the Rift Valley (ER) and West Africa (WA). Three models were then implemented as follows: The unconstrained null model (M0) allows geographic ranges to include any combination of areas (with input range and dispersal matrices in Lagrange left as the default values of 1.0). M1 restricts ranges to include a maximum of four adjacent areas, as follows: WA/EK/WK, ER/MTJ/EK, MTJ/EK/WK/SZ, EK/WK/SZ/SL, WK/SZ/SL, and SZ/SL, but does not restrict dispersal. The third more complex (stratified) model restricts ranges as in M1, and also constrains dispersal differentially across four time periods, 0-2, 2-5, 5-8 and 8-12 Mya, according to proposed changing patterns of drainage and physical barriers resulting from rifting. Thus, for the period 12-8 Mya, we constrain dispersal to occur only to and from regions SZ and SL, corresponding to the emergence of the common ancestor of Cryptomys/Fukomys in southern/South Africa (Faulkes et al., 2004;Ingram, Burda & Honeycutt, 2004). For the period 8-5 Mya, in our model dispersal is more relaxed as Fukomys spreads north (from SL and SK) through south central and east (WK, EK, MTJ and ER) to West Africa. Movement back again from WA to all regions, and from ER and MTJ to WA is not permitted during this phase (the latter reflecting proposed barriers formed by rifting); movement south from all regions to SL is prohibited except from SZ. As rifting and drainage evolution of major river systems (particularly the Zambezi and Kafue) is hypothesized to become more significant between 5-2 Mya (Van Daele et al., 2007a;Van Daele et al., 2007b;Moore, Cotterill & Eckardt, 2013), we further increase constraints on dispersal to also prohibit movement from SL and SZ north to WA, WK, EK, MTJ and ER. For the period 2-0 Mya and the final phases of rifting and drainage evolution, dispersal is only permitted from ER to MTJ, between MTJ and EK, from WK to SZ and SL, and between SZ and SL.

Morphology, craniometrics and analysis of skull shape
Pelage colour was recorded under natural daylight by consensus of three observers, with reference to Munsell Soil Color Charts (1954 Edition; Munsell Color Co., Inc. Baltimore, USA). Subsequent descriptions of colour all refer to this scale. Morphometric measurements were taken from a total of 34 skulls (Hanang/Mbulu region, n = 26; Ujiji, n = 5; Fukomys whytei, n = 3) using digital callipers (to the nearest 0.1 mm), as described by Verheyen, Colyn & Hulselmans (1996) (Fig. S1), together with standard head, body, tail and hind foot length measurements, in specimens where the entire body was available. Age classes were estimated from tooth eruption and wear characteristics as previously described for Fukomys damarensis (Bennett, Jarvis & Wallace, 1990).
Differences in craniometric measurements were investigated with MANOVA. Any measurements that were shown to have significant differences were then tested using a separate univariate ANOVA with a Tukey HSD to ascertain which species were statistically different. Bonferroni corrections were applied to these results to account for multiple testing. Because there were only three skulls available from F. whytei for craniometric measurements, this small sample size could potentially bias the statistical analysis. Therefore, a second MANOVA was carried out with just F. hanangensis and F. livingstoni measurements. A third MANOVA on just F. hanangensis craniometric measurements (the largest sample size) enabled us to look for sex and age class related differences in skull morphology. All ANOVA-based analyses were carried out using the 'stats' package in R 3.2.3 (R Core Team, 2015).
In order to investigate and quantify any differences in skull morphometrics between Fukomys livingstoni and Fukomys hanangensis, shape variation was investigated using landmark analysis, as previously described in Faulkes et al. (2010). The dorsal and ventral surfaces of skulls from specimens collected at Hanang (n = 2), Mbulu (n = 24), Ujiji (n = 5), and Kigogo (F. whytei, the geographically closest species; n = 3), were photographed three times to minimize the effects of misalignment. For further comparison, material (n = 20) from F. whytei was obtained in the form of photographs of the dorsal and ventral surfaces of the skulls collected by Allen & Loveridge (1933), and were provided by the Harvard Museum of Comparative Zoology (see Faulkes et al., 2010 for further information). In addition, photographs of the dorsal and ventral surfaces of the skulls from the more geographically remote F. anselli (n = 20) from Lusaka, Zambia were also included (this species is part of the micklemi clade in the genetic analysis). The relative locations of these samples are displayed in Fig. 1. To capture the shape, the 2-D coordinates of a total of 15 dorsal and 17 ventral landmarks as previously described in Faulkes et al. (2010;Fig. S2) were placed on each photograph and digitized using the TpsDIG2 software (Version 1.4; Rohlf, 2004), and mean relative warp scores for each specimen produced by tpsRelW (Version 1.36; Rohlf, 2003) were plotted.

New zoological taxonomic names
The electronic version of this article in Portable Document Format (PDF) will represent a published work according to the International Commission on Zoological Nomenclature (ICZN), and hence the new names contained in the electronic version are effectively published under that Code from the electronic edition alone. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank LSIDs (Life Science Identifiers) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix http://zoobank.org/. The LSID for this publication is: [urn:lsid:zoobank.org:pub:DC6D5104-CB60-48A1-9A06-B16A25DC6573]. The online version of this work is archived and available from the following digital repositories: PeerJ, PubMed Central and CLOCKSS.

General capture information
Ujiji A total of six animals were sampled from two sites 1.5 km apart and at an altitude of 2,601 to 2,624 m above sea level, at Msimba Village on the outskirts of Ujiji (Table 1; Fig. 1). The capture sites were either in or directly adjacent to fields containing maize, sweet potato, cassava, palms and bananas. At one location at Site 1 (Msimba village, Kasaka hamlet), an adult male and young female were captured from the same trap (NHMUK 2015.42 andNHMUK 2015.43), with an adult male (NHMUK 2015.46) trapped a few metres away and likely from the same burrow. An adult female and a young male (NHMUK 2015.44 andNHMUK 2015.45) were caught at further locations nearby (115 m distant) at Site 1 and may represent different colonies. All individuals had the same cyt b Haplotype ''D''. At Site 2 (Msimba village, Kabemba site), a single adult male (NHMUK 2015.47) was caught at 2,601 m above sea level in a valley with sweet potato fields and uncultivated grassland. This individual was assigned to a different cyt b Haplotype ''E''.

Hanang/Mbulu
At Hanang, 640 km east and slightly north of Ujiji, nine animals were collected from two sites around Hanang Mountain at altitudes of 1,856 to 1,957 m above sea level: at Gitting Village and at Jorodom Village, 8 km SW from Gitting (Table 1; Fig. 1). The capture sites were either in uncultivated grassland adjacent to or on the edge of fields next to cultivated mixed crops, bananas and planted trees. A single individual each had Haplotype ''A'' or Haplotype ''B'' at Gitting, while the remaining six at Gitting and two at Jorodom Village had Haplotype ''B'' (Table 1; Fig. 1).
Further north thirty-one animals were collected at altitudes of 2,115 to 2,188 m above sea level from seven distinct catching sites/colonies from another population at Mbulu, 42 km north of Gitting Village (Hanang). Habitat at Mbulu was similar to that around Hanang, either in uncultivated grassland adjacent to or on the edge of fields next to cultivated mixed crops (e.g., potatoes, maize, sugarcane, and bananas). From one to eight individuals were caught at the seven colony sites (Table 1). All twenty of the thirty-one animals sequenced from the Mbulu population had the same Haplotype ''C''.

Phylogenetic relationships
The maximum likelihood tree with the highest log likelihood (−5146.172) is shown in Fig. 2A. Initial trees for the heuristic search were obtained automatically by applying the maximum parsimony method. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.453)). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 52.275% sites; Tamura & Nei, 1993;Tamura et al., 2013). Bootstrap support was high for the main taxonomic groupings (80-100%), with the exception of Lufubu (69%), although the latter was consistently placed as the sister group to the East Bangweulu clade in all trees.
Maximum parsimony analysis produced four trees of length 883 ( Fig. 2B; consistency index = 0.466; retention index = 0.612). From the total of 1,140 characters, 674 were constant, 325 variable characters were parsimony informative and 141 uninformative. One of the trees (Fig. 2B (i)) was identical to the maximum likelihood tree. The other three most parsimonious trees differed in the relative placement of the whytei clade with respect to Lufubu, East Bangweulu, West Bangweulu and amatus clades, and the swapping of lineages of F. whytei from Mbala and Kigogo within the whytei clade ( Fig. 2B (ii), (iii) and (iv)).
Bayesian phylogenetic analysis performed with BEAST produced a tree identical to maximum likelihood tree 2a and maximum parsimony tree 2b (i), with all nodes having high posterior support (0.86 to 1.00/86-100%; Fig. 3).  Table 1, other species are designated according to current taxonomic understanding and GenBank Accession Numbers (Van Daele et al., 2007a;Van Daele et al., 2013;Faulkes et al., 2010).

Figure 3
Maximum clade credibility tree inferred using BEAST with an uncorrelated relaxed molecular clock model (allowing a variable rate of sequence evolution across the tree). The clock calibration to convert genetic distance to time is based on calibration 1 (Ingram, Burda & Honeycutt, 2004). The tree is also annotated with maximum likelihood reconstruction of geographic range evolution under the constrained (stratified) DEC model implemented in Lagrange (see text). Additional geological annotations below the tree are based on Voje et al. (2009), Trauth et al. (2005 and Ebinger et al. (1993), respectively. MTJ = corresponds to the period of increased tectonic activity at the Mbeya Triple Junction. Circled letters and the filled red circle on clade labels correspond to locations on Fig. 1, with GenBank accession numbers in parentheses after species names and locations. Geographic ranges of these extant taxa are denoted with abbreviations and colour coded boxes, and are defined below. Blue bars spanning nodes correspond to ages for the lower and upper bound of the 95% highest posterior density (HPD) intervals. Numbers at nodes refer to the posterior probabilities in support of that node (1 = 100%). Node A corresponds to the divergence of F. livingstoni, and Node B the divergence of F. hanangensis. The temporal position of these nodes, together with their HPD intervals, for clarity are further represented by the vertical dotted lines and the blue bars below the tree and above the geological annotation bar respectively The colour filled squares at nodes represent ancestral ranges for diverging lineages reconstructed by Lagrange. (continued on next page. . . )

Figure 3 (...continued)
The matrices above the tree associated with four time ranges depict possible area-to-area dispersals, with rows corresponding to ''from'' and columns ''to'', green fill ''yes'' and no fill ''no''. Areas as follows: SZ, south of the Zambezi River; SL, south of the Limpopo River; WK, west of the Kafue River; EK, east of the Kafue River; MTJ, Mbeya Triple Junction; ER, east of the Rift Valley; WA, West Africa. Relative positions of these geographic areas are depicted on the map inset, and follow the same abbreviations and colour codes as the filled squares on the tree nodes and taxon labels. Blue lines depict the Zambezi (Z), Kafue (K) and Limpopo Rivers (L), blue patches lakes, and bold black lines the main arms of the Rift Valley. Note that for clarity, other rivers are omitted.
Maximum likelihood, maximum parsimony and Bayesian trees were all congruent in supporting previously accepted taxonomic groupings (whytei, Lufubu, East Bangweulu, West Bangweulu, F. amatus, F. darlingi, F. damarensis, F. micklemi, F. bocagei, F. mechowii, and F. vandewoestijneae;Van Daele et al., 2007a;Van Daele et al., 2007b;Faulkes et al., 2010). The sequence of F. zechi from Ghana constitutes an extant representative of an early diverging lineage in the Fukomys clade, supporting previous studies that extra-limital populations of Fukomys in West and central Africa represent relic populations from an initial radiation of ancestral Fukomys (Ingram, Burda & Honeycutt, 2004). Next a clade containing F. mechowii, F. bocagei and the recently described F. vandewoestijineae (Van Daele et al., 2013) constitute a group of species distributed through central and west central Africa (Zambia and Angola), with an earlier common ancestor to the populations found at the geographically distant Ujiji. All trees consistently placed the animals collected from Ujiji and Hanang/Mbulu in reciprocally monophyletic groups, separated by the divergence of a major clade containing F. darlingi, F. damarensis, and F. micklemi/F. kafuensis. Finally, a monophyletic group containing five distinct clades was consistently recovered (West Bangweulu, F. amatus, Lufubu, East Bangweulu and F. whytei; Fig. 2).
The phylogenetic analysis provides evidence for two hitherto unrecognised phylogenetic species in Tanzania, one from Ujiji (Fukomys livingstoni sp. nov.) and a second from the Hanang/Mbulu region (Fukomys hanangensis sp. nov.). Both were genetically divergent from one another within the molecular phylogeny for the Fukomys genus, and also from the geographically closest clade containing F. whytei. Each of the Fukomys sp. nov. comprise clades strongly supported by bootstrap values of 100% (maximum likelihood) and posterior probabilities of 1.00 (Bayesian trees).

Inter-clade sequence divergence
Maximum Likelihood fits of 24 different nucleotide substitution models indicated that the Tamura-Nei + G + I (TN93+G+I) model of sequence evolution was the most appropriate. Both mean uncorrected-p and TN93+G+I corrected genetic distances between lineages/clades represented in Figs. 2 and 3 are displayed in Table 2. Uncorrected p (and (TN93+G+I) distances between F. hanangensis and ingroup lineages from different locations ranged from a minimum of 6.3% (7.0%) versus the Lufubu clade to 14.5% (18.5) versus F. zechi. Uncorrected p distances between F. livingstoni and ingroup lineages ranged from a minimum of 8.5% (9.8) versus F. amatus to 14.5% (18.8) versus F. zechi. Mean p distance between F. livingstoni and F. hanangensis was also 8.5% (9.9), while p distances between these and the geographically closest F. whytei clade were 9.7% (11.6) and 6.5%  Figure 3 summarises the molecular clock-based divergence times, together with 95% highest posterior density (HPD) intervals (equivalent to 95% confidence intervals), for the main nodes within the phylogeny generated using BEAST. According to the maximum clade credibility tree produced by BEAST (Fig. 3) Fig. 3), which forms the conduit between south central Africa and Tanzania. A sister group to F. hanangensis contains five clades with taxa restricted to Zambia, with the exception of the F. whytei clade (red circles in Fig. 1), which includes lineages that diverged much more recently in the Pleistocene. Geographically these populations of F. whytei are concentrated around the MTJ region, with only Kigogo population significantly within Tanzania. Assuming a confidence window of two log-likelihood units (Ree & Smith, 2008), Lagrange analysis of historical biogeography showed successively more likely scenarios when transitioning from the null (M0) to the more complex M1 and stratified models. This strongly supported our phylogeographic hypothesis, with the latter having a greater likelihood score (less negative) than the alternative M0 and M1 models (by 6.91 and 3.94 log-likelihood units respectively; Table 3). From the null model, estimates of the rate of dispersal successively increase while estimates of the rate of extinction decline. Thus, in the stratified model, the range data are better explained by a model of dispersal than a model of extinction. Figure 3 incorporates the results from Lagrange for the stratified model, and shows the optimal maximum likelihood reconstruction of geographic range evolution. This supports the hypothesis that from the common ancestor of Cryptomys and Fukomys, there was a dispersal north and range expansion into south-central and West Africa in Fukomys, with Cryptomys restricted to areas south of the Zambezi River (SZ). The lineages that follow the divergence of F. zechi into West Africa (WA) then show a pattern of cladogenesis with ancestral ranges for the descendant lineages encompassing Tanzania, east of the Rift Valley Table 3 Summary of inferences from Lagrange analysis of ancestral area and range evolution, under DEC models of historic biogeography. M0 represents the unconstrained (null) model that allows geographic ranges to include any combination of areas. M1 restricts ranges to include a maximum of three adjacent areas. The stratified model constrains dispersal differentially across four time periods, according to changing patterns of drainage and physical barriers resulting from rifting (Fig. 3).

Model
Log likelihood Fig. 4A (dorsal surface) and Fig. 4B (ventral surface). The clear separation of the F. livingstoni samples from the consensus shape (corresponding to the origin in the plot in Fig. 4A) and the F. hanangensis/F. whytei morphospace, along relative warps 1 and 2 of the dorsal surface is due to three main effects: (i) anterior-medial shifting of the jugal within the zygomatic arch (landmarks 4 and 6), (ii) shortening of the nasal bones, particularly at their posterior extent (landmarks 1, 2, 14 and 15), and (iii) anterior/anterior-medial shifts in the parietal (landmarks 8, 9 and 12; Fig. S2 and Fig. 4A). To some extent these same changes occur when moving from the F. hanangensis/F. whytei morphospace (i.e., upper right quadrant of Fig. 4A) to F. anselli (upper left quadrant of Fig. 4A), as this is also a shift along the x-axis. The additional separation of F. livingstoni from F. anselli points along the y-axis result from a posterior-lateral shifts in the narrowest inflection of squamosal (landmark 10) and the right anteriolateral tip of the parietal bone (landmark 9), and a posterior-medial shift in the anterior tip of the interparietal bone (landmark 11). A single small (35 g) male F. hanangensis is an outlier to the main cluster of points for this species, perhaps not being fully developed cranially.

Morphometric analysis of skulls
On the ventral skull surface, changes in F. livingstoni from the consensus shape (and F. hanangensis) are principally due to a small lateral shift in the posterior tip of auditory bulla at the junction with the occipital (landmark 9) and a small anterior-medial shift in the junction of the jugal and zygomatic process (landmark 4). These changes are present, Figure 4 Scatter plot of sample means of relative warp 1 (x-axis) against relative warp 2 (y-axis) from the shape analysis, together with thin plate splines showing landmark displacements from the consensus dorsal skull surface (the origin) to the extent of the diagonal plotted in each quadrant (e.g., upper left −0.14, 0.08 in (A)) for (A) dorsal, and (B) ventral skull surfaces. Plot symbols squares = males, circles = females, diamonds = sex unknown, F. hanangensis animals from Hanang. Colours: orange, F. livingstoni; black, F. hanangensis; grey, F. whytei; dark purple, F. anselli. Circled F. hanangensis outlier is animal 4332, a very small 35 g male (Table S1). but less exaggerated, in F. whytei, which also occupies the morphospace in the upper right quadrant of the relative warp plot (Fig. 4B). Changes from the consensus shape in the main group of F. hanangensis points were small as they cluster quite close to the origin on the plot. The separation of F. anselli along the y-axis results from anterior shifts in both the junction of squamosal and auditory bulla, and the lateral tip of auditory bulla (landmarks 7 and 8), and a posterior shift in the junction of jugal and zygomatic process (landmark 4; Fig. 4B).
Intraspecific differences in shape between the sexes was not apparent on either the dorsal or ventral plots for F. whytei, F. hanangensis or F. anselli, with points for males and females intermingled. With only a single female skull for F. livingstoni it is not possible to draw conclusions, although on both dorsal and ventral plots the female separates from the main group of male samples.
The full dataset of standard craniometric measurements of skulls are displayed in Table S1. Despite the small sample sizes for F. whytei and F. livingstoni, MANOVA analysis together with F. hanangensis revealed an overall global significant difference among the three species (Pillai's Trace = 1.84, F 2,28 = 3.59, P = 0.01), and significant differences in thirteen of the twenty-three measurements (summarized in Table S2). However, individual ANOVA on each of these significant measurements, followed by post-hoc Tukey tests with Bonferroni correction to investigate among species differences, left only one measurement that was significantly different (P < 0.05), due to the effects of the correction (for multiple testing) and small sample size for F. whytei. This metric was the greatest breadth of first upper molar (M13), with F. whytei larger than both F. hanangensis and F. livingstoni (2.30 ± 0.07 vs 1.90 ± 0.02 and 1.87 ± 0.04 mm respectively).
Because of the small sample size (n = 5) of F. livingstoni a comprehensive statistical analysis of craniometrics among sexes and age classes was not possible for this species. A MANOVA analysis of F. hanangensis including sex and age class as predictors revealed no significant differences in skull craniometrics between the sexes (Pillai's Trace = 0.80, F 1,19 = 15.04, P = 0.87) and among age classes (Pillai's Trace = 0.93, F 1,19 = 1.43, P = 0.54). The lack of sex differences in craniometric measurements is also evident in the shape analysis plots in Fig. 4, where the male and female points for F. hanangensis are intermingled.

Description of species
Family Bathyergidae Waterhouse, 1841 Genus Fukomys Kock et al., 2006 Fukomys livingstoni sp. nov. Livingstone's mole-rat (common name) LSID urn:lsid:zoobank.org:act:67DEACE5-3163-4FAE-885B-8EC04F072EEC Holotype NHMUK 2015.42 is an adult male collected in July 2013 at the Kasaka hamlet within the village of Msimba, near Ujiji. The specimen is composed of a skin and skull in very good condition (Figs. 5A and 5B; Fig. 6A). The external measurements (mm) are: head and body length 115.4, tail 8.9 and hind foot 22.0 ( Table 1). The body weight was 50 g. The pelage is darkish grey brown overall with a shorter very dark grey under-fur and a small irregularly shaped light grey head patch. Paratypes A further five specimens (NHMUK 2015.43 -NHMUK 2015 were collected from around the type locality (Table 1; Fig. 1). Etymology This species is named after Dr. David Livingstone, as Ujiji (the type locality) is the site of the famous meeting on 10 November 1871 when Henry Morton Stanley found the explorer David Livingstone, who many thought to be dead, and uttered the famous words ''Dr. Livingstone, I presume?'' (Stanley, 1872). Type locality Msimba village, 6.4 km northeast from the city centre of Ujiji (S04 • 51.760 ; E029 • 42.326 ). The specimen was trapped in a valley at an altitude of 793 m (2,601 ft) above sea level, in an area with moist sandy soil, where cassava, sweet potato, maize, palms and bananas were being cultivated.

Distribution and biology
The full range of this species remains to be determined with collection of the series described here restricted to around the village of Msimba on the outskirts of Ujiji. The holotype was captured from the same hole in the burrow as a young adult female (NHMUK 2015.43), with an adult male (NHMUK 2015.46) trapped a just few metres away and probably from the same burrow. The presence of adults and young adults together in a burrow suggests natal philopatry and cooperative breeding that is a characteristic of species within the genus Fukomys.

Diagnosis
Individuals of this species may be clearly distinguished from adjacent populations of F. hanangensis and F. whytei on the basis of morphology and molecular (DNA sequence) data. Morphologically, F. livingstoni is smaller with a mean adult body size (age class 2 and above) of 55 ± 8.9 g (n = 4) compared with F. hanangensis (mean adult body weight = 83.4 ± 5.6 g; n = 30; Table 1). Compared with F. hanangensis, the skull of F. livingstoni is shorter and narrower at the zygomatic arch, with a rostrum that is shorter and with less height at the mediosagittal projection at the anterior border of the first upper molars, but not different in breadth (Table S1 ; Fig. 6). A head spot (bles) is present in the specimens obtained for this study, although the presence/absence of the bles may not always be a reliable diagnostic feature as there is often specific variation in other bathyergids.

Description (and comparison with other species)
This is a small species of Fukomys: the four adults (age classes 2 to 4) ranged from 38-80 g in body weight of (mean = 55.0 ± 8.9 g; Table 1), similar in proportion to F. darlingi found in Zimbabwe, where mean body weights are 65.3 ± 14.1 g (males) and 62.9 ± 14.9 (females; Bennett, Jarvis & Cotterill, 1994). Sexual dimorphism in F. livingstoni remains to be investigated fully as there was only a single young female skull in our sample (age class 1; Table S1). When compared with other species of Fukomys from south-central Africa, the ratio of body length to body weight and the size of the skull (expressed as greatest width/greatest length) is at the lowest end of the distribution, and much smaller than F. damarensis and F. mechowii ( Fig. S4; Kingdon et al., 2013). Geometric morphometric analysis of the dorsal skull surface revealed differences between F. livingstoni and F. anselli, F. hanangensis and F. whytei. This is manifest as an anterior-medial shifting of the jugal within the zygomatic arch, shortening of the nasal bones, and anterior/anterior-medial shifts in the parietal. This results in the shorter narrower skull described in the diagnosis above, although not at the expense of rostrum width. On the ventral skull surface, changes observed in F. livingstoni were less pronounced, and principally due to a small lateral shift in the posterior tip of auditory bulla at the junction with the occipital and a small anterior-medial shift in the junction of the jugal and zygomatic process. Overall pelage colouration varied from darkish grey-brown to brown and dark brown, with shorter under-fur of very dark grey or black. A small irregularly-shaped head spot was present, varying from light grey to white in colour (Table S2; Fig. 5A). In this respect F. livingstoni is similar to neighbouring F. whytei, whose range extends into south-western Tanzania, where a small head spot is reported to be present in some populations (Allen & Loveridge, 1933;Burda et al., 2005). Otherwise F. livingstoni is much smaller than F. whytei, where body weight means (g) for animals from the type locality of Karonga, Malawi were 132.7 ± 22.3 (males, n = 4) and 121.5 ± 10.7 (females, n = 4; Burda et al., 2005). Specimens of F. whytei geographically closer to Ujiji (from Kigogo, Tanzania) were also within this size range, with a young animal of age class 1 (and 3/4 cheek teeth erupted) recorded at 101 g, larger than adult F. livingstoni (Table 1). To the trained eye, pelage colour may also distinguish the more grey-brown/brown F. livingstoni from F. whytei, although the latter is reportedly variable among populations from cinnamon-grey to dark slatey (Allen & Loveridge, 1933) and grey-buff (Burda et al., 2005). Body size also clearly distinguishes F. livingstoni from the larger F. hanangensis (see below) found further north in Tanzania, and while the latter lacks a paler coloured head spot (at least in the sample set reported here) and tends to be more yellowish brown, there is some overlap in pelage colouration.
Family Bathyergidae Waterhouse, 1841Genus Fukomys Kock et al., 2006 Fukomys hanangensis sp. nov. The Hanang mole-rat (common name) LSID urn:lsid:zoobank.org:act:59C00958-9628-461F-987D-AB897F52598F  (Table 1). The body weight was 62 g. The pelage is brown overall with a shorter black under-fur. No head spot is present. Paratypes A further 39 specimens including 27 paratypes (NHMUK 2015-14 and NHMUK 2015.16 -NHMUK 2015) and 12 samples retained at Queen Mary, University of London for further analysis. Eight of these were collected from around the type locality at Hanang, while the remaining 31 were from locations around Mbulu (Table 1; Fig. 1). Etymology Named after the location where the specimens were first collected around Mount Hanang, Tanzania. Type locality Jorodom village, (S04 • 29.510 ; E035 • 24.519 ). The specimen was trapped in a valley at an altitude of 1,957 m (6,422 ft) above sea level, in an uncultivated grass field surrounded by crop fields.

Distribution and biology
Currently the range of this species is known to be around Mount Hanang including the villages of Gitting and Jorodom, and extending to at least 40 km further north to Mbulu. The full range of this species remains to be determined. Aside from the first three animal captures in 2006, the remaining 37 specimens collected in 2009 from Gitting, Jorodom and Mbulu were gathered from 10 colonies with up to eight being caught from one burrow (Colony 4 at Tumati-Eyasirong, Mbulu; Table 1). These are not maximum colony sizes as burrows were not completely trapped out, and no breeding females were captured at Mbulu. Specimens from Colony 4 consisted of five males and three females, including a young 35 g male of age class 1, and mature adults of age classes 2 and 3. A similar spread of age classes was also seen among the animals collected from Mbulu Colony 5. These observations suggest natal philopatry and cooperative breeding for this species. Diagnosis Individuals of this species may be clearly distinguished from adjacent populations of F. livingstoni and the more geographically distant F. whytei on the basis of morphology and molecular (DNA sequence) data. Morphologically, F. hanangensis is larger than neighbouring F. livingstoni, with a mean adult body size (i.e., excluding animals of known age class 1) of 83.4 ± 5.6 g (range: 35-140; n = 30; Table 1). Compared with F. livingstoni, the skull of F. hanangensis is longer and wider at the zygomatic arch, with a rostrum that is longer and higher, but not broader (Table S1; Fig. 6).

Description (and comparison with other species)
F. hanangensis is a small to medium sized example of the genus Fukomys, while at an average adult size of 83 g it is larger than F. livingstoni (mean adult body weight = 55.0 ± 8.9 g, range: 38-80 g; n = 4; Table 1), it is slightly smaller in proportions to F. whytei, where body weight means (g) for animals from the type locality of Karonga, Malawi were 132.7 ± 22.3 (males, n = 4) and 121.5 ± 10.7 (females, n = 4; Burda et al., 2005). In comparison with other species of Fukomys from south-central Africa (Fig. S4), the ratio of body length to body weight and the size of the skull are smaller than F. damarensis and F. mechowii. However, it would be hard to distinguish F. hanangensis from species such as F. darlingi, F. anselli, F. bocagei, F. kafuensis and F. vandewoestijneae on the basis of body size alone. There was a trend towards male F. hanangensis being larger than females: 90.8 ± 8.0 g (n = 19) versus 72.5 ± 4.8 g (n = 10) respectively, although this was not significant (P = 0.127, t = 1.575, df = 27). However, there is no sexual dimorphism evident in either skull shape or craniometrics. Geometric morphometric analysis of the dorsal skull surface revealed no shape differences between F. hanangensis and F. whytei. However, when compared with F. anselli and F. livingstoni, there is a posterior-lateral shifting of the jugal within the zygomatic arch, lengthening of the nasal bones, and posterior/posterior-lateral shifts in the parietal. On the ventral skull surface, differences between F. hanangensis and F. livingstoni, F. whytei and (to a lesser extent) F. anselli, are more pronounced. This is manifest as a posterior shift in the posterior border of the molars, a small posterior-lateral shift in the junction of the jugal and zygomatic process, and a change in shape of the auditory bulla. The latter arises from posterior lateral shifts in both the junction of squamosal and auditory bulla, and the lateral tip of auditory bulla, and a small medial shift in the posterior tip of auditory bulla at the junction with the occipital. Overall pelage colouration varied from yellowish brown through dark yellowish brown to brown/dark brown. Under-fur was normally black, with two specimens very dark grey. None of the series described here had a lighter-coloured head spot present. Thus F. hanangensis can be morphologically distinguished from neighbouring F. livingstoni by its larger size and a more yellowy brown coat than F. livingstoni, and a lack of lighter head spot.

DISCUSSION
This study builds on preliminary sequence data from two mole-rats collected at Hanang in Tanzania, originally reported by Faulkes et al. (2010). Not only do we confirm the presence of a previously unrecognised species of African mole-rat in this region (Fukomys hanangensis), but also provide robust evidence for a second new species from specimens collected from a new locale at Ujiji (Fukomys livingstoni). We base our descriptions of the new species primarily on genetic data, although clear morphological differences are also evident in geometric morphometric analysis of skull shape, craniometrics, pelage and body size. Within the genus Fukomys, both F. hanangensis and F. livingstoni are distinct evolutionary lineages, as defined by the Phylogenetic Species Concept (Cracraft, 1989), and form separate clades nested among other clearly defined species in the cyt b molecular phylogeny. An increasing number of single locus mtDNA (in particular cyt b, but also 12S rRNA) phylogenies have produced robust and consistent phylogenies for the Bathyergidae and clarified the taxonomy of cryptic species (Allard & Honeycutt, 1992;Faulkes et al., 1997;Faulkes et al., 2004;Faulkes et al., 2010;Faulkes et al., 2011;Walton, Nedbal & Honeycutt, 2000;Huchon & Douzery, 2001;Ingram, Burda & Honeycutt, 2004). Recent phylogenomic analysis of the main lineages (genera) within the family using data from 3,999 concatenated genes agree fully with single gene studies, producing a tree with a congruent topology (Davies et al., 2015). Single nuclear gene-based studies (e.g., intron 1 of the nuclear transthyretin (TTR) gene) have also been useful, but such loci have been shown to be far less variable that mitochondrial cyt b, and thus lack phylogenetic signal and resolving power (see Ingram, Burda & Honeycutt, 2004). While Ingram et al. found that their TTR tree was not significantly different to their mitochondrial-based tree, the latter had more phylogenetic structure and the TTR data was unable to resolve most of the sub-clades within Fukomys. Importantly, it could not resolve such clearly defined species as F. damarensis, F. darlingi and F. kafuensis/micklemi, and therefore we did not sequence nuclear genes for this study. distinct radiations, perhaps also involving some local extinctions and replacements. In particular, F. livingstoni appears to be a lineage that dispersed into East Africa in the Pliocene (5.33-2.63 Mya), within a large clade with a common ancestor at Node A in Fig. 3 (the 95% HPD values for Node A are 4.89-2.58 Mya, according to the maximum clade credibility tree produced by BEAST). F. mechowi, F. bocagei and F. vandewoestijneae form an immediate outgroup to F. livingstoni, all of which are now distributed further south and west in Zambia, and into Angola in the case of F. bocagei (Fig. 1). F. hanangensis represents either a later, second incursion into Tanzania from a common ancestor occurring 3.25-1.68 Mya (Node B in Fig. 3), or the descendant from that common ancestor which also dispersed into Zambia. While Lakes Tanganyika and Nyasa may have formed a deep water barrier before the divergence of the F. hanangensis lineage and around the time of the earliest estimate for F. livingstoni, the dispersal route from south central Africa to East Africa was possible in the terrestrial corridor between Lakes Rukwa and Nyasa, that now forms the MTJ and Rungwe volcanic province (Branchu et al., 2005;Macgregor, 2015; see Fig. 1). The MTJ/Rungwe volcanic region is at the intersection of the Livingstone basin that forms the north east extremity of Lake Nyasa, the Rukwa-Songwe basin at the south east extremity of the Rukwa rift and the Usangu rift basin. While rifting and volcanism started in this area about 8.6 Mya and intra-basinal faulting, uplift and volcanism were particularly important in shaping the geology of the region at approximately 2.5 Mya, much of the modern topography was generated from 2 Mya and still continues to the present (Ebinger et al., 1993;Delvaux et al., 1998;Mortimer et al., 2007;Macgregor, 2015). Importantly, the range of timings for divergence of the common ancestors of both F. livingstoni and F. hanangensis (4.89-2.63 and 3.25-1.68 Mya respectively; Fig. 3) thus may precede the commencement of increased tectonic activity from 2.5 Mya to the present at the MTJ. It therefore seems highly likely that this increased faulting and uplift contributed to the separation of the south central populations of Fukomys and East African F. hanangensis and F. livingstoni, and today this mountainous habitat represents a significant physical barrier to dispersal for a subterranean rodent, with several points (e.g., Rungwe mountain) exceeding 2,900 m above sea level. Interestingly, F. whytei populations that diverged a little later are focussed mainly around the MTJ, with the exception of the Kigogo population having an earlier common ancestor that is slightly further east in Tanzania, perhaps arriving when dispersal was easier.
The phylogeographic analysis of F. hanangensis and newly acquired F. livingstoni samples add further support for our earlier assertions (Faulkes et al., 2010) that tectonic activity, climatic fluctuations and subsequent expansion and contraction of forest during the Pliocene-Pleistocene may have also played a role in the sub-structuring of populations and cladogenesis in Fukomys. The accompanying Pliocene expansion of C4 grasslands and the savannah habitat in this part of Africa, favoured by mole-rats, would likely have further facilitated range expansion of ancestral populations. The apparently localised and limited distribution of F. hanangensis and F. livingstoni in Tanzania makes assessment of their conservation status and other aspects of their biology a priority.
The Tanzania Commission for Science and Technology (COSTECH) granted a research permit for collection of rodents (permit no. 2013(permit no. -260-NA-2014 to Dr. Georgies Mgode.

DNA Deposition
The following information was supplied regarding the deposition of DNA sequences: The cytochrome b sequences are accessible via GenBank accession numbers KX905166 to KX905198.

Data Availability
The following information was supplied regarding data availability: The raw data has been supplied as Table S1.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.3214#supplemental-information.